"Harmonic oscillator 4" N = 4 a1 = zero(N,N) a2 = zero(N,N) for(k,0,N-1,for(j,0,N-1, a1[k+1,j+1] = sqrt(j + 1) (k == j + 1), a2[k+1,j+1] = sqrt(j) (k == j - 1) )) X = sqrt(hbar / (2 m omega)) (a1 + a2) P = i sqrt(m hbar omega / 2) (a1 - a2) "Commutator" dot(X,P) - dot(P,X) "Hamiltonian" dot(P,P) / (2 m) + 1/2 m omega^2 dot(X,X) C = 1/sqrt(2) (1,1,0,0) "Verify equation (1)" check(dot(conj(C),X,C) == sqrt(hbar / (2 m omega))) "ok" "Verify by integral method" psi(n) = C(n) 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) Psi = (psi(0) + psi(1)) / sqrt(2) xbar = integral(conj(Psi) x Psi, x) EXP = exp(-m omega x^2 / hbar) ERF = erf(sqrt(m omega / hbar) x) -- evaluate limits (minus infinity to plus infinity) xbar = eval(xbar, EXP, 0, ERF, 1) - eval(xbar, EXP, 0, ERF, -1) xbar check(xbar == sqrt(hbar / (2 m omega))) "ok" Run