Solution (source code)

= Solution

For a <bound state> take $b>0$ and set $\psi(r)=e^{-br}v(r)$. Substitution into the radial equation and cancellation of the common exponential give
$$
r v''+(2-2br)v'+(a-2b)v=0.
$$
Seek a <power series> $v(r)=\sum_{j\geq0}c_jr^j$ with $c_0\ne0$. Equating coefficients of $r^j$ yields the <polynomial construction of hydrogen S states> recurrence
$$
\boxed{c_{j+1}=\frac{2b(j+1)-a}{(j+1)(j+2)}c_j,\qquad j\geq0.}
$$
For any <integer> $N\geq1$, choose $b=a/(2N)$. The numerator vanishes at $j=N-1$, while it is nonzero for $j<N-1$, so $v$ is a nonzero <polynomial> of degree $N-1$. Therefore $\psi$ is continuous at $r=0$ and decays exponentially as $r\to\infty$. The three-dimensional normalization <integral> is
$$
4\pi\int_0^\infty r^2|\psi(r)|^2dr.
$$
It is finite near zero because $\psi$ is bounded, and finite at infinity because an exponential dominates every <polynomial>. Multiplying by its positive reciprocal square root gives a normalizable, nonzero <wavefunction>. The recurrence ensures the <differential equation> for $r>0$; at the singular endpoint the regular continuous extension has $2\psi'(0)+a\psi(0)=0$, so the potentially singular coefficient cancels. No nonexistent value of $1/r$ at zero is substituted.

Since $E=-\hbar^2b^2/(2m)$, the choice of $b$ gives
$$
\boxed{E_N=-\frac{\hbar^2a^2}{8mN^2}=-\frac{me^4}{32\pi^2\epsilon_0^2\hbar^2N^2},\qquad N=1,2,\ldots.}
$$
For $N=1$, the normalized example is $\psi(r)=\sqrt{b^3/\pi}\,e^{-br}$. For $N=2$, $v(r)$ is proportional to $1-br$. These illustrate the polynomial-times-exponential <Coulomb bound states> produced for every required energy. The radial derivative at the origin is generally nonzero, a Coulomb cusp consistent with the requested <continuity>.