  Magnetized neutron stars constitute a special class of compact objects
harbouring gravitational fields that deviate strongly from the Newtonian weak
field limit. Moreover strong electromagnetic fields anchored into the star give
rise to non-linear corrections to Maxwell equations described by quantum
electrodynamics (QED). Electromagnetic fields close to or above the critical
value of $\BQ=4.4\times10^9$~T are probably present in some pulsars and for
most of the magnetars. To account properly for emission emanating from the
neutron star surface like for instance thermal radiation and its polarization
properties, it is important to include general relativistic (GR) effects
simultaneously with non-linear electrodynamics. This can be achieved through a
3+1 formalism known in general relativity and that incorporates QED
perturbations to Maxwell equations. Starting from the lowest order corrections
to the Lagrangian for the electromagnetic field, as given for instance by
Born-Infeld or Euler-Heisenberg theory, we derive the non-linear Maxwell
equations in general relativity including quantum vacuum effects. We also
derive a prescription for the force-free limit and show that these equations
can be solved with classical finite volume methods for hyperbolic conservation
laws. It is therefore straightforward to include general relativity and quantum
electrodynamics in the description of neutron star magnetospheres by using
standard classical numerical techniques borrowed from Maxwell and Newton
theory. As an application, we show that spin-down luminosity corrections
associated to QED effects are negligible with respect to GR corrections.
