"Hydrogen atom 2"

H(f) = -hbar^2 D(f) / (2 mu) - hbar^2 / (mu a0 r) f

D(f) = 1/r^2 d(r^2 d(f,r), r) +
       1/(r^2 sin(theta)) d(sin(theta) d(f,theta), theta) +
       1/(r sin(theta))^2 d(f,phi,2)

psi(n,l,m) = R(n,l) Y(l,m)

R(n,l) = 2 / n^2 *
         a0^(-3/2) *
         sqrt((n - l - 1)! / (n + l)!) *
         (2 r / (n a0))^l *
         L(2 r / (n a0), n - l - 1, 2 l + 1) *
         exp(-r / (n a0))

-- associated Laguerre polynomial (k is a local var)

L(x,n,m,k) = (n + m)! sum(k, 0, n, (-x)^k / ((n - k)! (m + k)! k!))

-- spherical harmonic

Y(l,m) = (-1)^m sqrt((2 l + 1) / (4 pi) (l - m)! / (l + m)!) *
         P(l,m) exp(i m phi)

-- associated Legendre of cos theta (arxiv.org/abs/1805.12125)

P(l,m,k) = test(m < 0, (-1)^m (l + m)! / (l - m)! P(l,-m),
           (sin(theta)/2)^m sum(k, 0, l - m,
           (-1)^k (l + m + k)! / (l - m - k)! / (m + k)! / k! *
           ((1 - cos(theta)) / 2)^k))

E(n) = -1 / (2 mu) (hbar / a0)^2 / n^2

E1 = E(1)
E2 = E(2)
E3 = E(3)

"Verify eigenfunctions"

-- n = 1, 2, 3, ...
-- l = 0, 1, ..., n - 1
-- m = -l, ..., 0, ..., l

f = psi(1,0,0)
check(H(f) == E1 f)

f = psi(2,0,0)
check(H(f) == E2 f)

f = psi(2,1,0)
check(H(f) == E2 f)

f = psi(2,1,1)
check(H(f) == E2 f)

f = psi(2,1,-1)
check(H(f) == E2 f)

"ok"