Solution (source code)

= Solution

For a <barotropic fluid>, $\nabla p=\rho\nabla h$. Magnetostatic force balance and $L\alpha=F(\beta)$ give
$$
\rho\nabla(\Phi+h)=\mathbf F_m
=\frac{L\beta-FF'}{\mu_0r^2}\nabla\beta.
$$
The left side is a gradient multiplied by $\rho$, so taking the curl shows that $(L\beta-FF')/(r^2\rho)$ is constant on each $\beta$ surface. Absorbing the fixed factor $\mu_0$ into an arbitrary function $G$ gives the <Axisymmetric magnetostatic Grad-Shafranov system>
$$
\boxed{L\beta=F(\beta)\frac{dF}{d\beta}+r^2\rho G(\beta)}.
$$