Solution (source code)

= Solution

For bound <hydrogen atom> states, the radial <wave function> is regular and finite at the origin, and decays at infinity so that $4\pi\int_0^\infty r^2|\Psi(r)|^2\,dr<\infty$. Equivalently the reduced radial function $w(r)=r\Psi(r)$ has $w(0)=0$ and is square integrable on the half-line. A singular $1/r$ branch is not an admissible regular Coulomb <eigenfunction>; its <derivative> energy diverges and it introduces an origin singularity absent from the <Hamiltonian> domain.

Write $\kappa=e^2/(4\pi\epsilon_0)$, so $a=\hbar^2/(m\kappa)$. Multiplication by $r$ turns the radial equation into
$$
w''+\left(-\mu^2+\frac2{ar}\right)w=0,qquad
\mu^2=-2mE/\hbar^2.
$$
Since $w=a e^{-\lambda R}f(R)$ and $r=aR$, substitution gives
$$
\boxed{f''-2\lambda f'+(\lambda^2-\mu^2a^2+2/R)f=0.}
$$
Choose $\boxed{\lambda=\mu a>0}$ to remove the constant coefficient and factor out the decaying exponential. For $f=\sum_{k\geq0}a_kR^{k+c}$, the lowest power has coefficient $c(c-1)a_0$, so $\boxed{c(c-1)=0}$. The $c=0$ behavior would give $\Psi\sim a_0/R$, violating regularity; choose $c=1$.

Matching successive powers then yields
$$
\boxed{a_{k+1}=\frac{2\{\lambda(k+1)-1\}}{(k+2)(k+1)}a_k.}
$$
A nonterminating tail eventually has one sign, and its coefficient ratio is asymptotic to $2\lambda/k$. More precisely, for any $0<\delta<\lambda$, the ratio is eventually at least $(2\lambda-\delta)/(k+1)$. Its same-sign tail consequently dominates a positive constant times the exponential series of rate $2\lambda-\delta$, up to finitely many terms. The factor $e^{-\lambda R}$ cannot then make $\Psi$ decay. Thus square integrability forces termination.

If the last nonzero term is $a_{n-1}R^n$, termination requires $\lambda n=1$. Therefore
$$
\boxed{E_n=-\frac{\hbar^2}{2ma^2n^2}=-\frac A{2n^2},\qquad
A=\frac{\hbar^2}{ma^2}=\frac{me^4}{16\pi^2\epsilon_0^2\hbar^2}.}
$$
For $n=1$, $\lambda=1$, $f$ is proportional to $R$, and
$$
\boxed{\Psi_0(r)\propto e^{-r/a}.}
$$
The radial volume element is $4\pi r^2dr$, so the mean radius is the ratio of the radius-weighted <norm> <integral> to the <norm> <integral>. With $r=aR$, the common angular and <normalization> factors cancel:
$$
\langle r\rangle=a\frac{\int_0^\infty R^3e^{-2\lambda R}\,dR}
{\int_0^\infty R^2e^{-2\lambda R}\,dR}
=\frac{3a}{2\lambda}.
$$
In the <ground state> $\boxed{\langle r\rangle=3a/2}$.