A numerical code using Hermite bicubic finite elements has been developed for the computation of axisymmetric magnetohydrodynamic (MHD) equilibria. The code provides a mapping to flux coordinates for MHD stability calculations. Several lest cases are studied to show the convergence rate for the equilibrium. Convergence tests are also presented for the eigenvalues of the stability calculations when the equilibrium mesh is varied.