Solution (source code)

= Solution

Put $I=RB_\phi$ and define $\Delta_*\psi=\psi_{RR}-R^{-1}\psi_R+\psi_{zz}$. Direct cylindrical differentiation gives
$$
(\nabla\times B)_R=-I_z/R,\quad
(\nabla\times B)_z=I_R/R,\quad
(\nabla\times B)_\phi=-\Delta_*\psi/R.
$$
The azimuthal component of the <Lorentz force> is proportional to $I_z\psi_R-I_R\psi_z$. A <force-free magnetic field> therefore has parallel meridional gradients of $I$ and $\psi$, implying $I=f(\psi)$ on connected regular flux surfaces. Consequently
$$
\boxed{B_\phi=f(\psi)/R.}
$$
With this dependence the poloidal part of $\nabla\times B$ is $f'(\psi)B_p$. Its full cross product with $B$ is
$$
(\nabla\times B)\times B
=-\frac{\Delta_*\psi+f(\psi)f'(\psi)}{R^2}\nabla\psi.
$$
Vanishing force therefore requires $\Delta_*\psi+ff'=0$ wherever the flux gradient is nonzero, extending by regularity through ordinary isolated nulls. Finally, evaluating the three-dimensional cylindrical divergence gives $R^2\nabla\cdot(R^{-2}\nabla\psi)=\Delta_*\psi$, so
$$
\boxed{R^2\nabla\cdot(R^{-2}\nabla\psi)+f\frac{df}{d\psi}=0.}
$$
This derives the <Grad-Shafranov equation for a force-free magnetic field>. The function $f$ is set by boundary/current information; it is not an additional fixed constant. Its single-valued global form presumes the usual connected flux-surface labeling. The equation also makes $\nabla\times B=f'(\psi)B$, explicitly showing that current is parallel to the <magnetic field>.