Soft Gamma-Ray Repeaters and Anomalous X-Ray Pulsars are extreme manifestations of the most magnetized neutron stars: magnetars. The phenomenology of their emission and spectral properties strongly support the idea that the magnetospheres of these astrophysical objects are tightly twisted in the vicinity of the star. Previous studies on equilibrium configurations have so far focused on either the internal or the external magnetic field configuration, without considering a real coupling between the two fields. Here we investigate numerical equilibrium models of magnetized neutron stars endowed with a confined twisted magnetosphere, solving the general relativistic Grad-Shafranov equation both in the interior and in the exterior of the compact object. A comprehensive study of the parameters space is provided to investigate the effects of different current distributions on the overall magnetic field structure.