 Solution

The flux-function form solves $\nabla \cdot \mathbf{B}=0$ identically. With $k=k(\psi )$, the poloidal relation $\rho {\mathbf{u}}_p=k{\mathbf{B}}_p$ and axisymmetry solve steady mass conservation. The ideal induction equation is solved by the velocity representation with the field-line angular velocity $\omega =\omega (\psi )$. The azimuthal momentum equation integrates once to the magnetohydrodynamic angular-momentum invariant $\ell =\ell (\psi )$.

What remains is the pressure or entropy equation and the two poloidal components of momentum. Projecting poloidal momentum along a field line produces the Bernoulli invariant in part (c); projecting across magnetic surfaces produces the transfield or Grad-Shafranov equation that determines their shape. An equation of state and boundary conditions complete the problem. The stellar gravitational potential is prescribed, so no self-gravity Poisson equation remains to solve.