"Harmonic oscillator 1"

H(f) = -hbar^2 / (2 m) d(f,x,x) + 1/2 m omega^2 x^2 f

psi(n) = C(n) exp(-i E(n) t / hbar) *
         exp(m omega x^2 / (2 hbar)) *
         d(exp(-m omega x^2 / hbar), x, n)

C(n) = (-1)^n / sqrt(2^n n!) *
       (m omega / (pi hbar))^(1/4) *
       (hbar / (m omega))^(n/2)

E(n) = hbar omega (n + 1/2)

psi0 = psi(0)
psi1 = psi(1)
psi2 = psi(2)
psi3 = psi(3)
psi4 = psi(4)

E0 = E(0)
E1 = E(1)
E2 = E(2)
E3 = E(3)
E4 = E(4)

"Verify wavefunctions"

check(i hbar d(psi0,t) == H(psi0))
check(i hbar d(psi1,t) == H(psi1))
check(i hbar d(psi2,t) == H(psi2))
check(i hbar d(psi3,t) == H(psi3))
check(i hbar d(psi4,t) == H(psi4))

check(E0 psi0 == H(psi0))
check(E1 psi1 == H(psi1))
check(E2 psi2 == H(psi2))
check(E3 psi3 == H(psi3))
check(E4 psi4 == H(psi4))

-- integral of |f|^2 from minus infinity to plus infinity

I(f) = 2 eval(integral(conj(f) f, x),
       exp(-m omega x^2 / hbar), 0,
       erf(sqrt(m omega / hbar) x), 1)

check(I(psi0) == 1)
check(I(psi1) == 1)
check(I(psi2) == 1)
check(I(psi3) == 1)
check(I(psi4) == 1)

-- time-independent wavefunctions

psi0 = psi0 / exp(-i E0 t / hbar)
psi1 = psi1 / exp(-i E1 t / hbar)
psi2 = psi2 / exp(-i E2 t / hbar)
psi3 = psi3 / exp(-i E3 t / hbar)
psi4 = psi4 / exp(-i E4 t / hbar)

check(E0 psi0 == H(psi0))
check(E1 psi1 == H(psi1))
check(E2 psi2 == H(psi2))
check(E3 psi3 == H(psi3))
check(E4 psi4 == H(psi4))

"ok"