  Using an axisymmetric numerical code, we perform an extensive study of the
magnetic field configurations in non-rotating neutron stars, varying the mass,
magnetic field strength and the equation of state. We find that the monopolar
(spherically symmetric) part of the norm of the magnetic field can be described
by a single profile, that we fit by a simple eighth-order polynomial, as a
function of the star's radius. This new generic profile applies remarkably well
to all magnetized neutron star configurations built on hadronic equations of
state. We then apply this profile to build magnetized neutron stars in
spherical symmetry, using a modified Tolman-Oppenheimer-Volkov (TOV) system of
equations. This new formalism produces slightly better results in terms of
mass-radius diagrams than previous attempts to add magnetic terms to these
equations. However, we show that such approaches are less accurate than usual,
non-magnetized TOV models, and that consistent models must depart from
spherical symmetry. Thus, our `universal' magnetic field profile is intended to
serve as a tool for nuclear physicists to obtain estimates of magnetic field
inside neutron stars, as a function of radial depth, in order to deduce its
influence on composition and related properties. It possesses the advantage of
being based on magnetic field distributions from realistic self-consistent
computations, which are solutions of Maxwell's equations.
