Abstract: Accurate tokamak plasma equilibrium solution in flux coordinates is crucial for many stability and transport studies. Different approaches for dealing with singularities in solving the nonlinear Grad–Shafranov (GS) equation in flux coordinates or also known as straight field line coordinates are proposed in this paper. The GS equation is solved by iterating the position of grids directly in flux coordinates, and hence, no additional errors are introduced due to mapping process for a convergent solution. The singularity at magnetic axis in flux coordinates is removed by using a novel coordinate transform technique. Different from other techniques previously developed, no assumption in boundary condition at magnetic axis is used. This is consistent with the fact that there is no physical boundary at the magnetic axis. A flux coordinate system with poloidal coordinate chosen as the geometric poloidal angle is proposed. It conquers the difficulty in no definition of poloidal coordinate in flux coordinates at separatrix because of the singularity at x-point(s) in a divertor configuration. It also simplifies the process for computing poloidal flux coordinate during the iteration for solving the nonlinear GS equation. Non-uniform grids can be applied in both radial and poloidal coordinates, which allows it to increase the spacial resolution near x-point(s) in a divertor configuration. Based on the model proposed in this paper, a new Flux coordinates based EQuilibrium solver (FEQ) in tokamaks is developed. The numerical solutions from this code agree well with both the analytic Solov’ev solution and the numerical one from the EFIT code for a divertor configuration in the EAST tokamak. This code can be applied for simulating different equilibria with prescribed shape, pressure and current profiles, i.e. including both limiter and divertor configurations, positive triangularity and negative triangularity, different β, arbitrary magnetic shear profile etc. It provides a powerful and convenient fixed-boundary inverse equilibrium solver including both magnetic axis and separatrix in the solution for tokamak researches.