A method is presented for the calculation of an axisymmetric toroidal magnetohydrodynamic equilibrium in a flux coordinate. The algorithm is based on a combination of a useful finite element method (FEM) and an iteration scheme. The mesh structure for calculation is corrected for the finite elements always to be along magnetic surfaces. High accuracy is thus attained in the vicinity of the magnetic axis and in the regions where magnetic surfaces are closely located. Results produced by this code compared with a special solution. The method is useful for obtaining equilibria to use in stability and transport calculation.