(Part of the Wolverhampton Lectures of Physics's Quantum Physics Course)
However accurate, Bohr's model was however clearly incomplete and handwaving. In particular, although it provides correctly the energy spectrum of Hydrogen, this was for this particular case only, for which it works exceedingly well, but even there, it finds quickly its limitation, for instance, it does not say anything about the details of the transitions between the level (the rate at which they occur and the underlying physical mechanism).
The full picture is given by solving Schrödinger's equation for the wavefunction of the electron:
$$H\psi=E\psi$$
where the Hamiltonian is the electron's kinetic energy in the Coulomb potential
$$H=-{\hbar^2\over 2m}\nabla^2-{e^2\over4\pi\epsilon_0}{1\over r}$$
Note that doing so, we do not describe the proton's dynamics, which only provides the potential. This assumes that the proton has an infinite mass as compared to the electron, which can be relaxed by working in the reference frame not of the proton itself but of the center of mass of the compound system, in which case, $m\rightarrow mM/(m+M)$ where $M$ is the mass of the proton.
Since $H$ is spherically symmetric, the problem is solved for the wavefunction $\psi(r,\theta,\phi)=R(r)Y(\theta,\phi)$, with $Y$ a spherical harmonic and the corresponding radial equation for $u(r)=rR(r)$ reads:
$$-{\hbar^2\over2m}{d^2u\over dr^2}+\left[-{e^2\over4\pi\epsilon_0}{1\over r}+{\hbar^2\over 2m}{l(l+1)\over r^2}\right]u=Eu\,.$$
Introducing the usual $\kappa$ number, this time for negative energies $E<0$ as the problem will have both bound and scattering states:
$$\kappa={\sqrt{-2mE}\over\hbar}$$
and also introducing a rescaled radius:
$$\rho\equiv\kappa r\quad\text{where}\quad\rho_0\equiv{me^2\over2\pi\epsilon_0\hbar^2\kappa}$$
so that the radial equation reads:
$${d^2u\over d\rho^2}=\left[1-{\rho_0\over\rho}+{l(l+1)\over\rho^2}\right]u\,.\tag{1}$$
To solve this equation, we consider first the asymptotics:
$$u(\rho)\rightarrow A\exp(-\rho)\,.$$
$${d^2u\over d\rho^2}={l(l+1)\over\rho^2}u\,.$$
$$u(\rho)\rightarrow C\rho^{l+1}\,.$$
We now introduce still another variable $v(\rho)$, to "peel off" the asymptotic behaviour:
$$u(\rho)=\rho^{l+1}e^{-\rho}v(\rho)\,.$$
If we now compute the successive derivatives, we find:
$$\begin{align} {du\over d\rho}&=\rho^le^{-\rho}\left[(l+1-\rho)v+\rho{du\over d\rho}\right]\,,\\ {d^2u\over d\rho^2}&=\rho^le^{-\rho}\left\{\left[-2l-2+\rho+{l(l+1)\over\rho}\right]v+2(l+1-\rho){dv\over d\rho}+\rho{d^2v\over d\rho^2}\right\}\,, \end{align}$$
so that, inserting this into the radial equation (1), we find the equation in terms of $v$:
$$\rho\frac{d^2v}{d\rho^2}+2(l+1-\rho)\frac{dv}{d\rho}+[\rho_0-2(l+1)]v=0\,.\tag{2}$$
The reason to bring it in this form is that we can now solve it with a series expansion for $v$:
$$v(\rho)=\sum_{j=0}^\infty a_j\rho^j$$
for which, also computing the respective derivatives:
$$\begin{align} \frac{dv}{d\rho}&=\sum_{j=0}^\infty ja_j\rho^{j-1}=\sum_{j=0}^\infty(j+1)a_{j+1}\rho^j\,,\\ {d^2v\over d\rho^2}&=\sum_{j=0}^\infty j(j+1)a_{j+1}\rho^{j-1}\,, \end{align}$$
which, once inserted into Eq. (2), yields
$$\sum_{j=0}^\infty\left\{j(j+1)a_{j+1}\rho^j+2(l+1)(j+1)a_{j+1}\rho^j-2ja_j\rho^j+[\rho_0-2(l+1)]a_j\rho^j\right\}=0\,.$$
Setting the coefficients with the same power equal to zero, which is the only way to cancel the sum, we find a recursive relation:
$$a_{j+1}=\frac{2(j+l+1)-\rho_0}{(j+1)(j+2l+2)}a_j\,.\tag{3}$$
Let us consider again the asymptotic behaviour. For $j$ large enough, Eq. (3) becomes $a_{j+1}\approx(2/j)a_j$ with solution, neglecting the departures for small $j$
$$a_j={2^j\over j!}a_0\tag{4}$$
in which case, the solution would be $v(\rho)=a_0\sum_{j=0}^\infty(2\rho)^j/j!=a_0e^{2\rho}$, but since $u(\rho)=A\rho^{l+1}e^{-\rho}v(\rho)$, we have again an exponential divergence as $\rho\rightarrow\infty$. At such, the solution is also not physical. Unless... the series is finite, in which case the asymptotic behaviour of Eq. (4) can be excluded. This therefore brings a new constrain:
$$2(j+l+1)-\rho_0=0\,.$$
The only variable here is $j$, so there exist, for each $l$, a $j_\mathrm{max}$ that cancels the numerator of the recursive relation, i.e., such that $j_\mathrm{max}+l+1=\rho_0$. This brings another integer, the principal quantum number $n$ (which makes the connection to the quantization introduced by Bohr):
$$n\equiv j_\mathrm{max}+l+1$$
at which point:
$$\rho_0=2n\,.$$
In this way, we have linked $n$ to the energy, since $\rho_0$ was itself energy dependent. This yields:
$$E_n=-{m\over2\hbar^2}\left({e^2\over4\pi\epsilon_0}\right)^2{1\over n^2}$$
and therefore:
$$E_n={E_1\over n^2}$$
where $E_1$ is as defined in the old Bohr theory, but it has now be derived as part of the Schrödinger equation. This also provides a sound derivation for the Bohr radius $a_\mathrm{B}$ which is linked to $\rho$ as
$$\rho={r\over a_\mathrm{B}n}\,.$$
We can now work out recursively the wavefunctions for the various quantum numbers by plugging back $v(\rho)$ in the radial wavefunction
$$R_{nl}(r)={1\over r}\rho^{l+1}e^{-\rho}v(\rho)\,.$$
for the various quantum numbers $n=1, 2, \ldots$ with $j_\mathrm{max}=n-l-1$ for all the possible positive integer values of $l$. In terms of $l$, since $l=n-1-j_\mathrm{max}$, the leading term of the series ($a_0$) always provides a solution on its own, when $l=n-1$, in which case $j_\mathrm{max}=0$ and the series is a single term. Since $l\ge0$, $j_\mathrm{max}$ is at most $n-1$, there are a maximum of $n$ terms in the series for the $l=0$ wavefunction. In between these two extremes, we have series $n$ values of $l$, which we can all work out in turns:
$${\displaystyle R_{10}(r)={a_0\over a_\mathrm{B}}\exp(-r/a_\mathrm{B})}\,.$$
$$\int_0^\infty|R_{10}|^2r^2\,dr={a_0^2\over a_\mathrm{B}^2}\int_0^\infty e^{-2r/a_\mathrm{B}}r^2\,dr={a_0^2a_\mathrm{B}\over4}=1$$
$$\boxed{\displaystyle\psi_{100}(r,\theta,\phi)={1\over\sqrt{\pi a_\mathrm{B}^3}}\exp(-r/a_\mathrm{B})\,.}$$
$$a_1=-a_0\quad\text{and}\quad a_2=0$$
$$R_{20}(r)={a_0\over2a_\mathrm{B}}\left(1-{r\over 2a_\mathrm{B}}\right)e^{-r/2a_\mathrm{B}}\,.$$
$$R_{21}(r)={a_0\over4a_\mathrm{B}^2}re^{-r/(2a_\mathrm{B})}\,.$$
For each $l$ there is a corresponding number of $m$ solutions (the third quantum number from spherical harmonics). In the general case, from $j=n-l-1$, we find that for each $n$, there are $n$ possible values of $l$, namely:
$$\boxed{l=0, 1, 2, \ldots, n-1}$$
And since for each $l$ we have $(2l+1)$ values of $m$, the degeneracy $d(l)$ (i.e., the number of solutions) of the energy level $E_n$ is
$$d(l)=\sum_{l=0}^{n-1}(2l+1)=n^2\,.$$
As for the series, it turns out that it can be explicited in closed-form in terms of the Laguerre polynomials $L_q$ defined as
$$L_q(x)\equiv e^x\left({d\over dx}\right)^q(e^{-x}x^q)\,.$$
From this, one can define the so-called associated Laguerre polynomials
$$L_{q-p}^p(x)\equiv(-1)^p\left({d\over dx}\right)^pL_q(x)\,,$$
in which case, the function $v(\rho)$ is given by:
$$v(\rho)=L^{2l+1}_{n-l-1}(2\rho)\,.$$
The first few radial wavefunctions can be worked out to be:
\begin{align} R_{10}(r)&=2a_\mathrm{B}^{-3/2}\exp(-r/a_\mathrm{B})\\ R_{20}(r)&={1\over\sqrt{2}}a_\mathrm{B}^{-3/2}\left(1-{1\over2}{r\over a_\mathrm{B}}\right)\exp(-r/(2a_\mathrm{B}))\\ R_{21}(r)&={1\over\sqrt{24}}a_\mathrm{B}^{-3/2}{r\over a_\mathrm{B}}\exp(-r/(2a_\mathrm{B})) \end{align}
The normalization can be found (we are not deriving it but give the result directly):
$$\boxed{\psi_{nlm}(r,\theta,\phi)=\sqrt{\left({2\over na_\mathrm{B}}\right)^3{(n-l-1)!\over2n(n+l)!^3}}\exp\left(-{r\over na_\mathrm{B}}\right)\left({2r\over na_\mathrm{B}}\right)^lL_{n-l-1}^{2l+1}\left({2r\over na_\mathrm{B}}\right)Y_l^m(\theta,\phi)\,.}$$
These functions thus defined form an orthonormal set.
The visualisation of these wavefunctions is a famous problem that we will tackle on the next lecture.