The code XNS solves for the axisymmetric equilibrium configuration of differentially rotating Neutron Stars with toroidal magnetic field, using the XCFC approximation for the metric, in spherical coordinates.  The code is based on the metric module and routines developed for the GR-MHD code X-ECHO

You can download an older version of the code  XNS here !
And an older Guide!

There are several other codes published in the years (RNS, RotStar, etc....) you might want to check.
However, solutions, even for rapid rotators close to the mass shedding limit, show that the difference
between conformally flat metric CFC/XCFC, and the more appropriate quasi-isotropic coordinates is
order of 0.1% (larger differences are expected for strongly distorted cases of disk-like configurations), so in principle one would expect the CFC/XCFC limit to provide a reasonably good approximation of the
correct solution for Neutron Stars. There are some computational benefits in using XCFC over quasi-isotropic coordinates, and  as an Astrophysicist I am looking more toward a solution that is physically acceptable than one that is mathematically correct.

These are some comparisons with the code RNS, for a uniform rotator close to mass shedding limit. You can see that the errors are ~0.1% or smaller. On the left is a comparison of the polar, and equatorial densities: the solid lines are from RNS and the dashed lines (can you see them?)  from XNS. On the right is a comparison of the errors in the metric terms. All the lines that you see represent various metric terms (see the Guide for explanation). They are all < 0.1%, the dash dotted line is a measure of how far from conformally-flat the RNS result is. That indicates the quality of the CFC/XCFC approximation.

Here is instead a differentially rotating magnetized solution.

There are still a few things that can/must be improved:
  • The code is limited to Politropic EoS. It should not be hard to implement more realistic and/or user supplied EoS. But this needs some adjustments in a few routines.
  • The code needs at the moment some hand-made search for the initial starting density (see the Guide for explanation). The search could be automated, but I have not yet thought how to do it cleverly, due to the oscillatory nature of the convergence. 
  • The metric solver utilizes routines developed for a more general axisymmetric geometry, without symmetry with respect to the equator. This also could be optimized for Neutron Stars.
  • The code does not handle poloidal magnetic fields.  However, equilibrium with poloidal fields is strictly only possible for uniform rotators, with poloidal magnetic field fully confined within the star. Unless the Neutron Star is in a true vacuum, which can never be realized/maintained in reality, the solution with a weak poloidal field, will at best be quasi-equilibrium. Adding a weak poloidal component is one of the future possible upgrades, but it will require some major changes which will not be coming soon.
 ....... need more here! coming soon!