%Chapt6Exercise4c.m
%simple harmonic oscillator wave function psi_n(xi)and 
%modulus wavefunction squared |psi_n(xi)|^2
%xi=(m*w/hbar)*x
%using relation 
%psi_n=((2/n)^0.5)*((xi*psi_n-1)-((((n-1)/2)^0.5)*psi_n-2))
clear;
clf;
npoints=500;			%(number of data points in plot) - 1
nlim=100;				%arbitrary limit to value of n that can be plotted

n=input('Input quantum index n = ');%Input from keyboard value of quantum index n
									%of wave function to be plotted (n=0,1,2,3, ... )

if n < 0;    error('minimum value of n must be greater or equal to 0');     end;
if n > nlim; error('maximum value of n must be less than or equal to 100'); end;

ximax=sqrt((2*n)+1);					%classical turning point
xiplot=(3/((n+1)^(1/3)))+(1.2*ximax);	%plot range of x-axis
deltaxi=2*xiplot/npoints;				%increment in xi

Ao=(1/pi)^(1/4);						%Normalization amplitude for n=0
An=1.1*Ao/((n+1)^0.1);					%fix vertical scale

xi(1)=-xiplot;							%first value of xi

for j=2:1:npoints+1
   xi(j)=xi(j-1)+deltaxi;
   psi1(j)=Ao*exp((-xi(j)^2)/2); 	%known n=0 ground-state wave function
   psi2(j)=(sqrt(2))*xi(j)*psi1(j);	%known n=1 first excited-state wave function
   psi(j)=psi2(j);
end   
	if n < 1; psi=psi1; end;
   
if n>=2
      for ni=2:1:n
        nin=double(ni);
for j=2:1:npoints+1
   psi(j)=(sqrt(2/nin))*((xi(j)*psi2(j))-((sqrt((nin-1)/2))*psi1(j)));
   psi1(j)=psi2(j);					%update new value of psi_(n-2)
   psi2(j)=psi(j);					%update new value of psi_(n-1)
end    
		end
end

figure(1); 
subplot(1,2,1),plot(xi,psi);
axis([-xiplot,xiplot,-An,An]);
xlabel('Position, \xi (m)'),ylabel('Wave function, \psi(\xi)');
ttl = sprintf('Chapt6Exercise4c, xi =(m*omega/hbar)*x, n=%3.0f',n);
title (ttl);
subplot(1,2,2),plot(xi,abs(psi.^2));
axis([-xiplot,xiplot,0,An^2]);
xlabel('Position, \xi (m)');
ylabel('Modulus wave function squared, |\psi(\xi)| ^2');
ttl2=sprintf('Classical turning-point=+/-%5.3f',ximax);
title (ttl2);
