  We derive general equations for axisymmetric Newtonian MHD and use these as
the basis of a code for calculating equilibrium configurations of rotating
magnetised neutron stars in a stationary state. We investigate the field
configurations that result from our formalism, which include purely poloidal,
purely toroidal and mixed fields. For the mixed-field formalism the toroidal
component appears to be bounded at less than 7%. We calculate distortions
induced both by magnetic fields and by rotation. From our non-linear work we
are able to look at the realm of validity of perturbative work: we find for our
results that perturbative-regime formulae for magnetic distortions agree to
within 10% of the nonlinear results if the ellipticity is less than 0.15 or the
average field strength is less than $10^{17}$ G. We also consider how
magnetised equilibrium structures vary for different polytropic indices.
