-- Verify Dirac's equation (32)

-- See "Quantum Mechanics and a Preliminary Investigation of the Hydrogen Atom"

psi(n,m) = 1 / sqrt(pi a0 (n + 1/2)) *
           sqrt((n - abs(m))! / (n + abs(m))!) *
           (2 r / a0 / (n + 1/2))^abs(m) *
           L(2 r / a0 / (n + 1/2), n - abs(m), 2 abs(m)) *
           exp(-r / a0 / (n + 1/2)) *
           exp(i m phi)

a0 = hbar^2 / (k mu)

-- Laguerre polynomial

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

-- energy eigenvalues

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

-- integrator

I(f) = do(
  f = defint(f, phi, 0, 2 pi),
  f = integral(f,r),
  0 - eval(f,r,0)
)

N = 5 -- arbitrary cutoff for demo

R = zero(N,N)
Enm = zero(N,N)

for(n,1,N, for(m,1,N, R[n,m] = I(conj(psi(n,0)) r psi(m,0))))

for(n,1,N, Enm[n,n] = E(n))

-- solve for P

P = sqrt(-2 mu / hbar^2 (Enm + k inv(R)))

print(R,P)

"Verify Dirac's equation (32)"

Pn = P + n hbar unit(5,5)

H(P) = -hbar^2 / (2 mu) inv(dot(P,P))

check(H(Pn) - H(P) == hbar^2 / (2 mu) (inv(dot(P,P)) - inv(dot(Pn,Pn))))

"ok"

"Verify Schroedinger equation"

check(-hbar^2 / (2 mu) dot(P,P) - k inv(R) == Enm)

"ok"