```%
%  Lecture 6  Harmonic oscillator
%

home
clear

m=1;
k=1;
hbar=1;
omega=(k/m)^0.5;

alpha=(m*k/hbar/hbar)^0.25;

norm=@(n) (alpha/(2^n*factorial(n)*pi^0.5))^0.5;

harmonic=@(n,x) norm(n)*hermite(n,alpha*x).*exp(-(alpha*x).^2/2);
E = @(n) (n+0.5)*hbar*omega;

Lim=3.5;
x=[-Lim:0.05:Lim];

h0 = harmonic(0,x);
h1 = harmonic(1,x);
h2 = harmonic(2,x);
h3 = harmonic(3,x);
h4 = harmonic(4,x);

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked');

plot(x,0.5*k*x.*x);
hold all

plot(x,h0+E(0),[-Lim Lim],[E(0) E(0)],'k:')
plot(x,h1+E(1),[-Lim Lim],[E(1) E(1)],'k:')
plot(x,h2+E(2),[-Lim Lim],[E(2) E(2)],'k:')
plot(x,h3+E(3),[-Lim Lim],[E(3) E(3)],'k:')
plot(x,h4+E(4),[-Lim Lim],[E(4) E(4)],'k:')
text(-3,E(0)+0.1,'n=0')
text(-3,E(1)+0.1,'n=1')
text(-3,E(2)+0.1,'n=2')
text(-3,E(3)+0.1,'n=3')
text(-3,E(4)+0.1,'n=4')

xlabel('distance'); ylabel('Energy')
title('Harmonic oscillator functions')
axis([-3.5 3.5 0 5.25])

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked');

plot(x,0.5*k*x.*x);
hold all

plot(x,h0.*h0+E(0),[-Lim Lim],[E(0) E(0)],'k:')
plot(x,h1.*h1+E(1),[-Lim Lim],[E(1) E(1)],'k:')
plot(x,h2.*h2+E(2),[-Lim Lim],[E(2) E(2)],'k:')
plot(x,h3.*h3+E(3),[-Lim Lim],[E(3) E(3)],'k:')
plot(x,h4.*h4+E(4),[-Lim Lim],[E(4) E(4)],'k:')
text(-3,E(0)+0.1,'n=0')
text(-3,E(1)+0.1,'n=1')
text(-3,E(2)+0.1,'n=2')
text(-3,E(3)+0.1,'n=3')
text(-3,E(4)+0.1,'n=4')

xlabel('distance'); ylabel('Energy')
title('Harmonic oscillator functions squared')
axis([-3.5 3.5 0 5.25])

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked');
hold all
x=[-5:0.05:5];
h7 = harmonic(7,x);
plot(x,h7.*h7,[-5 5],[0 0],'k:')
% classical turning point:
Xclassic=sqrt((2*7+1)*hbar/(k*m)^0.5);
plot([-Xclassic -Xclassic],[-0.6 0.6],'r:')
plot([Xclassic Xclassic],[-0.6 0.6],'r:')
axis([-5 5 0 0.3])
title('n=7 Harmonic oscillator probability, including classical turning points')
xlabel('Distance');ylabel('Probability (1/distance)');
```