-- Compute diagonal elements of H in equation (10) of

-- "Accurate energies of the He atom with undergraduate quantum mechanics"

-- by Robert C. Masse and Thad G. Walker

-- https://pages.physics.wisc.edu/~tgwalker/118.HeAtom.pdf

-- equation (4)

P(r,n,l) = sqrt(2 (n - l - 1)! / n^2 / (n + l)!) *
           (4 r / n)^(l + 1) *
           exp(-2 r / n) *
           L(4 r / n, n - l - 1, 2 l + 1)

-- Laguerre polynomial

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

-- equation 5

E0(n1,n2) = -2 / n1^2 - 2 / n2^2

-- R is distance
-- N is number of intervals

R = 20.0
N = 400
dr = R / N

f1 = zero(N)
f2 = zero(N)
f3 = zero(N)
f4 = zero(N)

for(k, 1, N, f1[k] = P((k - 0.5) dr, 1, 0)^2)
for(k, 1, N, f2[k] = P((k - 0.5) dr, 2, 0)^2)
for(k, 1, N, f3[k] = P((k - 0.5) dr, 3, 0)^2)
for(k, 1, N, f4[k] = P((k - 0.5) dr, 4, 0)^2)

-- numerical integrator

I(u,v) = sum(j, 1, N, sum(k, 1, N,
 u[j] v[k] / (test(j > k, j, k) - 0.5)
))

H11 = E0(1,1) + I(f1,f1) dr
H22 = E0(1,2) + I(f1,f2) dr
H33 = E0(1,3) + I(f1,f3) dr
H44 = E0(1,4) + I(f1,f4) dr

H11
H22
H33
H44

-- compare with equation (10)

check(abs(H11 - (-2.750)) < 0.001)
check(abs(H22 - (-2.080)) < 0.001)
check(abs(H33 - (-2.023)) < 0.001)
check(abs(H44 - (-2.010)) < 0.001)

"ok"