Solution (source code)

= Solution

Apply the <Frobenius method> at the origin. The lowest power in $f''-\ell(\ell+1)f/r^2$ is $r^{c-2}$, with coefficient
$$
\bigl[c(c-1)-\ell(\ell+1)\bigr]a_0.
$$
The terms $-2\alpha f'$ and $2mKZf/(\hbar^2r)$ start at power $r^{c-1}$ and therefore cannot cancel a nonzero $r^{c-2}$ coefficient. Because $a_0\ne0$, the <indicial equation> is
$$
\boxed{c(c-1)-\ell(\ell+1)=0,\qquad (c-\ell-1)(c+\ell)=0.}
$$
Its roots are $c=\ell+1$ and $c=-\ell$. Since $\psi=r^{-1}f e^{-\alpha r}$, their leading behaviors are respectively $\psi\sim r^\ell$ and $\psi\sim r^{-\ell-1}$. The physical choice is
$$
\boxed{c=\ell+1.}
$$
This is the <Coulomb radial origin condition>. For $\ell\geq1$, the other branch is not even locally normalizable: its radial <probability> integral behaves as $\int_0^\varepsilon r^2r^{-2\ell-2}\,dr=\int_0^\varepsilon r^{-2\ell}\,dr$, which diverges.

For $\ell=0$, the singular behavior $1/r$ is locally <square-integrable> with the measure $r^2dr$, so local <wavefunction normalization> alone would not reject it. The ordinary three-dimensional <Coulomb potential> problem additionally requires the regular <boundary condition> $u(0)=0$. A $1/r$ branch violates that condition and has divergent <kinetic energy> integral $\int_0^\varepsilon r^2|\partial_r(1/r)|^2\,dr=\int_0^\varepsilon r^{-2}\,dr$. Thus it is excluded from the physical energy domain as well. The regular choice gives $u\sim r^{\ell+1}$ in every case.