In [1]:

    

using PyPlot


    

In [2]:

    

function H(x, n)
    if n < 0
        error("negative n not allowed")
    elseif n == 0
        return Float64(1)
    elseif n == 1
        return 2*x
    else
        memory = Array{Float64}(3)
        memory[1] = 1
        memory[2] = 2*x
        for i in 2:n
            memory[3] = 2*x*memory[2]-2*(i-1)*memory[1]
            memory[1] = memory[2]
            memory[2] = memory[3]
        end
        return memory[2]
    end
end
function phi(x,n)
    return pi^(-1/4)/sqrt(2^n*factorial(n))*H(x,n)*exp(-1//2*x^2)
end
function V(x)
    return 1//2*x^2
end
function E(n)
    return n+1//2
end

function trapeze(func, start, stop, N, args=[])
    h = (stop-start)/(N-1)
    I = -h*(func(start,args...)+func(stop,args...))/2
    for x = linspace(start, stop, N)
        I += h*func(x, args...)
    end
    return I
end
;


    

In [3]:

    

figure(1)
X = Array(linspace(-5,5,201))
plot(X, V.(X), label="Potential")
for n in 0:4
    plot(X, phi.(X,n)+E(n), label="\$\\bar\\phi_$n\$")
end
xlabel("x")
ylabel("Energie")
title("Quantenmechanischer harmonischer Oszillator")
legend()
show()


    

In [4]:

    

function prod_phi_mn(x,m,n)
    return phi(x,m)*phi(x,n)
end
D = Array{Float64}(5,5)
for m in 0:4
    for n in 0:4
        D[m+1,n+1] = trapeze(prod_phi_mn,-5,5,100000,[m,n])
    end
end
D


    

    
Out[4]:





5×5 Array{Float64,2}:
  1.0           2.8563e-18   -5.54028e-11  -1.30139e-18  -7.51719e-10
  2.8563e-18    1.0          -2.80557e-17  -1.5994e-9     1.86331e-17
 -5.54028e-11  -2.80557e-17   1.0           1.38514e-17  -2.71767e-8 
 -1.30139e-18  -1.5994e-9     1.38514e-17   1.0          -1.30348e-17
 -7.51719e-10   1.86331e-17  -2.71767e-8   -1.30348e-17   1.0        

In [ ]: