% this is a good idea it erases
clear ;


% Create a grid.   
%
% I. The grid runs from xstart... xstop.
%
% II. dx is the grid spacing.
%
% III. The grid is 
%
%      xstore(1) = xstart, x_store(2) = xstart+dx, .... x_store(N) = xstop
%
% IV. N is the length of the grid
xstart  = -5.  ;
xstop   = 5.  ;
dx      = 0.01 ;
x_store = xstart:dx:xstop ;   
N       = length(x_store) ;

% Store the potential for plotting
k = 1. ;
hw0 = 1. ;
for i = 1:N
   x = x_store(i) ;
   v_store(i) = 0.5*k*x^2  ;
end

% Choose the energy
e = 0.5*hw0

% This  bit of code determines the location
% of the classical turning points
x1 = -sqrt(2.*e/k) ;
x2 = sqrt(2.*e/k) ;


%Plot the potential and the energy line. Also draw the 
% classical turning points
%
% The first command tells matlab to plot in the top window 
% the potential and the energy line.
% 
% You should change xmin1 xmax1 and ymax1 as you need.
subplot(2,1,1), plot(x_store,v_store) 
ymax1 = 10 ;
xmin1 = -4.  ;
xmax1 = 4.;
axis([xmin1 xmax1 0.0 ymax1 ] )

% Plot the energy line
line([xmin1, xmax1], [e, e], 'Linestyle', '-.') 

% These commands plot the turning point lines
line([x1,x1], [0. ymax1], 'LineStyle', '--')
line([x2,x2], [0. ymax1], 'LineStyle', '--')

% These commands add the title etc
title('Potential V(x),  \Delta\epsilon_{0} =\pm 0.003') 
ylabel({'V(x) /[(h/2\pi)^2 /M L^2]' })
xlabel({'x/L' })
str1 = sprintf('\\epsilon = %7.4f, \\Delta x= %5.3f', e, dx)  ;
text(0.7,0.9,str1,'Units', 'normalized') 
 

% create an array to store the wave function
psi_store = zeros(1,N) ;

% Integrate from  xmin up to zero
%
% We start at some small value at xmin  
psi  = 1.e-4 ;    % The wave function at the current point 
dpsi = 0. ;       % The derivative of the wave fucntion at the current point

psi_store(1) = psi ;
for i = 1:N-1
   % Determine where we are
   x = x_store(i)  ;
   % Update the wave by one grid spacing
   %
   % psi (x + dx) = psi  + dx*dspi
   % dpsi (x + dx) = psi + dx*ddpsi
   %
   % the second derivative ddpsi is given by the
   % Schrodinger equation
   psi = psi + dx*dpsi ; 
   v = 0.5*k*x*x ;
   ddpsi = -2.*(e - v) * psi ;
   dpsi = dpsi +  dx*ddpsi ;

   % Store the wave function at psi(x + dx) in the array for later plotting
   psi_store(i+1) = psi  ;
end


% Plot the wave function in the lower sub-pannel
subplot(2,1,2), plot(x_store, psi_store) 

% You should change this as you  need
% These parameters are the limits of the graph
ymin =-15. ;
ymax = 15.  ;
xmin = -4 ;
xmax = 4 ;
axis([xmin xmax ymin ymax] ) ;

% Print the tuning lines on the lower graph
line([x1,x1], [ymin ymax], 'LineStyle', '--')
line([x2,x2], [ymin ymax], 'LineStyle', '--')

% Print the title and other info
title('Wave function \psi_{0}(x),  \Delta\epsilon_{0} = \pm 0.003') 
xlabel('x/L')
ylabel('\psi(x) (arb units)')
str = sprintf('\\epsilon = %7.4f, \\Delta x= %5.3f', e, dx)  ;
text(0.7,0.9,str,'Units', 'normalized')