%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% harmonic_pot.m
%
% Matlab Scrödinger Equation solution a particle in harmonic potential.
%
% for simplicity we took hbar = m = 1
%
%
% - Dec 2016
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

V = inline('x*x/2'); % the potential
L = 1.0;       % width
n = 2000;      % number of parts
xL= -10*L;     % left  end where psi(xL)=0
xR= +10*L;     % right end where psi(xR)=0
dx= (xR-xL)/n;
x = xL;
A = zeros(n-1); % matrix elements

for i=1:n-1
   for j=1:n-1
     if (i==j+1 || i==j-1)
         A(i,j) = -0.5/dx^2;
     end
   end
   x = x + dx;
   A(i,i) = 1.0/dx^2 + V(x);
end

[wfn E] = eig(A); % get eigenvalues and eigenvectors

x = xL+dx:dx:xR-dx;
fprintf('*** Energy eigenvalues ***\n',n);

for i=1:4
   energy(i) = E(i,i);
   wfn(:,i) = wfn(:,i)/max(wfn(:,i)); % normalize the wfn to 1
   fprintf('n=%2d --> E=%15.4f\n',i,energy(i));
   plot(x,wfn(:,i).*wfn(:,i),'*-');
   axis([xL xR 0.0 3.0]);
   pause(1);
end
