% Plot the first few 1D harmonic oscillator wavefunctions

x=[-4:.01:4];
y = zeros(8,801);
par = x.^2/2;

for n=0:7
    y(n+1,:) = exp(-x.^2/2) .* polyval(HermitePoly(n),x) / (pi^(1/4)*sqrt(2^n*factorial(n)))/2 ...
    + n + 1/2;
end

figure;
plot(x,par,'k',...
    x,1/2,x,y(1,:),x,3/2,x,y(2,:),x,5/2,x,y(3,:),x,7/2,x,y(4,:),x,9/2,x,y(5,:),x,11/2,x,y(6,:),...
    x,13/2,x,y(7,:),x,15/2,x,y(8,:)),
    xlabel('x'),...
    ylabel('\psi_n(x)'),title('1D Harmonic Oscillator Wavefunctions (n\leq 7)')
grid on;