  Extremely magnetized neutron stars with magnetic fields as strong as $\sim
10^{15-16}$ G, or magnetars, have received considerable attention in the last
decade due to their identification as a plausible source for Soft Gamma
Repeaters and Anomalous X-ray Pulsars. Moreover, this class of compact objects
has been proposed as a possible engine capable of powering both Long and Short
Gamma-Ray Bursts, if the rotation period in their formation stage is short
enough (~1 ms). Such strong fields are expected to induce substantial
deformations of the star and thus to produce the emission of gravitational
waves. Here we investigate, by means of numerical modeling, axisymmetric static
equilibria of polytropic and strongly magnetized stars in full general
relativity, within the ideal magneto-hydrodynamic regime. The eXtended
Conformally Flat Condition (XCFC) for the metric is assumed, allowing us to
employ the techniques introduced for the X-ECHO code [Bucciantini & Del Zanna,
2011, Astron. Astrophys. 528, A101], proven to be accurate, efficient, and
stable. The updated XNS code for magnetized neutron star equilibria is made
publicly available for the community (see
www.arcetri.astro.it/science/ahead/XNS). Several sequences of models are here
retrieved, from the purely toroidal (resolving a controversy in the literature)
or poloidal cases, to the so-called twisted torus mixed configurations,
expected to be dynamically stable, which are solved for the first time in the
non-perturbative regime.
