-- psi returns a hydrogen atom eigenfunction psi(n,l,m) = R(n,l) Y(l,m) -- R returns a radial eigenfunction 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!)) -- Bohr radius a0 = 4 pi epsilon0 hbar^2 / (e^2 mu) -- Y returns a spherical harmonic eigenfunction Y(l,m) = (-1)^m sqrt((2l + 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), 2^(-m) sin(theta)^m sum(k,0,l - m, (-1)^k (l + m + k)! / (l - m - k)! / (m + k)! / k! * ((1 - cos(theta)) / 2)^k)) -- H is the Hamiltonian H(psi) = -hbar^2 Lap(psi) / (2 mu) - e^2 / (4 pi epsilon0 r) psi Lap(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) -- E(n) returns the nth energy eigenvalue E(n) = -mu / (2 n^2) (e^2 / (4 pi epsilon0 hbar))^2 "Verify eigenfunctions" check(H(psi(1,0,0)) == E(1) psi(1,0,0)) check(H(psi(2,0,0)) == E(2) psi(2,0,0)) check(H(psi(2,1,0)) == E(2) psi(2,1,0)) check(H(psi(2,1,1)) == E(2) psi(2,1,1)) check(H(psi(2,1,-1)) == E(2) psi(2,1,-1)) check(H(psi(3,0,0)) == E(3) psi(3,0,0)) check(H(psi(3,1,0)) == E(3) psi(3,1,0)) check(H(psi(3,1,1)) == E(3) psi(3,1,1)) check(H(psi(3,1,-1)) == E(3) psi(3,1,-1)) check(H(psi(3,2,0)) == E(3) psi(3,2,0)) check(H(psi(3,2,1)) == E(3) psi(3,2,1)) check(H(psi(3,2,-1)) == E(3) psi(3,2,-1)) check(H(psi(3,2,2)) == E(3) psi(3,2,2)) check(H(psi(3,2,-2)) == E(3) psi(3,2,-2)) "ok" Run