Solution (source code)

= Solution

In an axisymmetric magnetostatic state, pressure and gravitational forces are poloidal, so the azimuthal component of the <Lorentz force density> must vanish. Hence $f=0$ and
$$
\nabla(L\alpha)\times\nabla\beta=0.
$$
The two gradients are locally parallel, so $L\alpha$ is constant on each regular level surface of $\beta$. Therefore
$$
\boxed{L\alpha=F(\beta)}
$$
for an arbitrary flux function $F$.