%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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