Tokamak MHD equilibrium plays a crucial role in magnetic confinement nuclear fusion studies. It is the foundation of many hot research topics, such as stability, transport, heating, and current drive by radio frequency waves. Due to toroidal symmetry, the MHD equilibrium problem reduces to a 2D boundary value problem of nonlinear Grad-Shafranov (G-S) equation. Typically, solving the tokamak MHD equilibrium requires solving the G-S equation boundary value problem on an irregular domain, which brings difficulties to numerical solving. This paper presents a novel finite difference method-harmonic mesh finite difference method-for solving 2D parabolic boundary value problems on arbitrary irregular domains. This method consists of two parts: mesh generation and solution of problem. Since mesh generation can be regarded as transforming an irregular domain into a rectangular domain, and the transformation adopted is a harmonic transformation, the resulting mesh is called a harmonic mesh. On the harmonic mesh, the boundary problem is solved numerically by using the finite difference method. Since the problem is solved based on the harmonic mesh, numerical accuracy can be guaranteed. Based on the method, a simulation platform on tokamak MHD equilibrium has been developed, named TESP (Tokamak Equilibrium Simulation Platform), which includes four parts: input module, mesh generation module, equilibrium solving module, and output module. The platform has been validated by using the Solov'ev analytical solution, and it has been used to provide the spatial distributions of the equilibrium configurations. Finally, several MHD equilibrium configurations for CFETR (China Fusion Engineering Experimental Reactor) and ITER (International Thermonuclear Experimental Reactor) are obtained by using TESP. (c) 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license(https://creativecommons.org/licenses/by/4.0/).