The Grad-Shafranov equation describes the magnetic flux distribution of plasma in an axisymmetric system like a tokamak-type nuclear fusion device. The equation is transformed into an equivalent boundary integral equation by expanding the inhomogeneous term related to the plasma current into a polynomial. In the present research, the singularity of the fundamental solution, which consists of two elliptic integrals, and the properties of singular integrals have been minutely investigated. The discontinuous quadratic boundary elements have been introduced to give an accurate solution with a small number of boundary elements. Ampere's circuital law has been applied to estimate the total plasma current from the boundary integral of the poloidal field, and based on this idea, a new scheme to calculate the eigenvalue has also been proposed.3