A new type of boundary element method (BEM) has been applied to solve the Grad-Shafranov equation and to give a distribution of magnetic flux function in a tokamak nuclear fusion device. The quantity mu(0)rjphi related to the plasma current profile is expanded into a two-dimensional polynomial. Using the particular solution of the Grad-Shafranov equation with this inhomogeneous polynomial source and applying Green's second identity, the domain integral related to the plasma current is transformed into an equivalent boundary integral. Domain discretization is not required in this formulation, thus preserving all the advantages of the BEM. Numerical computations of all boundary integrals are only required in the initial stage of the eigenvalue iteration, so that the number of eigenvalue iterations does not hamper the total computing time. Test calculations demonstrated that the present method provides stable and accurate solutions.