  Magnetized neutron stars are privileged places where strong electromagnetic
fields as high as $\BQ=4.4\times10^9$~T exist, giving rise to non-linear
corrections to Maxwell equations described by quantum electrodynamics (QED).
These corrections need to be included to the general relativistic (GR)
description of a magnetic dipole supposed to be anchored in the neutron star.
In this paper, these QED and GR perturbations to the standard flat space-time
dipole are calculated to the lowest order in the fine structure
constant~$\alpha_{\rm sf}$ and to any order in the ratio $\Rs/R$ where $R$ is
the neutron star radius and $\Rs$ its Schwarzschild radius. Following our new
3+1~formalism developed in a previous work, we compute the multipolar
non-linear corrections to this dipole and demonstrate the presence of a small
dipolar~$\ell=1$ and hexapolar~$\ell=3$ component.
