Solution (source code)

= Solution

\b[The fundamental nonnegative mode is $\phi_0J_0(2.405r/R)\cos(\pi z/H)$.] Put the origin at the centre of the <cylinder>, so $-H/2\leq z\leq H/2$. Neglecting extrapolation imposes zero <neutron flux> at $r=R$ and $z=\pm H/2$. Write
$$
B_m^2=\frac{(\eta-1)\Sigma_a}{D},
\qquad
\phi(r,z)=F(r)Z(z).
$$
For an axisymmetric field, <separation of variables> in <cylindrical coordinates> gives
$$
\frac{F''+r^{-1}F'}F+\frac{Z''}Z+B_m^2=0.
$$
Choose the axial separation constant so $Z''+\beta^2Z=0$. Midplane symmetry and the lowest positive axial mode give $Z(z)=\cos(\beta z)$, with
$$
\beta H/2=\pi/2,\qquad \beta=\pi/H.
$$
The radial equation is
$$
F''+\frac1rF'+\alpha^2F=0,
\qquad \alpha^2=B_m^2-\beta^2.
$$
The solution regular on the axis is the zero-order <Bessel function of the first kind>, $F(r)\propto J_0(\alpha r)$; the independent second-kind solution is singular there. The first radial zero gives $\alpha R=j_{0,1}\simeq2.405$. Higher zeros would introduce sign changes and do not describe the fundamental nonnegative <neutron flux>. Normalize by the central value, using $J_0(0)=1$:
$$
\boxed{\phi(r,z)=\phi_0J_0(\alpha r)\cos(\beta z),
\qquad \alpha\simeq\frac{2.405}{R},\quad\beta=\frac{\pi}{H}.}
$$
The printed $2.405$ is a rounded <Bessel function> zero, rather than an exact root.