Solution (source code)

= Solution

Put
$$
q(R)=\frac{4\pi\alpha(R)}c.
$$
For an <axisymmetric vector field> depending only on the <cylindrical radius> $R$, the <solenoidal vector field> condition is
$$
\frac1R\frac{d}{dR}(RB_R)=0.
$$
Hence $RB_R$ is constant. The hypothesis $B_R\to0$ on the axis forces
$$
\boxed{B_R=0}.
$$
The azimuthal and axial components of $\nabla\times\mathbf B=q\mathbf B$ are then
$$
-\frac{dB_z}{dR}=qB_\phi,
\qquad
\frac1R\frac{d}{dR}(RB_\phi)=qB_z.
$$
Eliminating $B_\phi=-B_z'/q$ gives the closed <ordinary differential equation>
$$
\boxed{B_z''+\left(\frac1R-\frac{q'}q\right)B_z'+q^2B_z=0},
\qquad q=\frac{4\pi\alpha}{c}.
$$
Once $B_z$ is known, $B_\phi=-B_z'/q$ and $B_R=0$ determine the full <cylindrical force-free magnetic field>.