Solution (source code)

= Solution

The remaining <magnetic field> has $\mathbf B=B_\phi(R)\mathbf e_\phi+B_z(R)\mathbf e_z$. The radial <Lorentz force density> separates into a <magnetic pressure> gradient and inward <magnetic tension>:
$$
\frac1{4\pi}[(\nabla\times\mathbf B)\times\mathbf B]_R
=-\frac1{8\pi}\frac{d(B_z^2+B_\phi^2)}{dR}-\frac{B_\phi^2}{4\pi R}.
$$
Here $(\mathbf B\cdot\nabla)\mathbf B$ has radial component $-B_\phi^2/R$ because $\partial_\phi\mathbf e_\phi=-\mathbf e_R$. Thus <magnetostatic equilibrium> gives
$$
\frac d{dR}\left(p+\frac{B_z^2+B_\phi^2}{8\pi}\right)+\frac{B_\phi^2}{4\pi R}=0.
$$
Integrating the axial component of the magnetostatic <Ampère-Maxwell equation> gives $I(R)=cRB_\phi/2$, since regularity removes the integration constant. Substituting $B_\phi=2I/(cR)$ combines the toroidal terms into
$$
\boxed{\frac d{dR}\left(p+\frac{B_z^2}{8\pi}\right)
+\frac1{2\pi c^2R^2}\frac{dI^2}{dR}=0.}
$$
This is the <cylindrical magnetostatic pressure balance> equation in terms of enclosed <electric current>.