% -1/2 d2psi/dxi2 + V psi = epsilon psi
% calculations for matrix A:
% =-1/2 (fj+1 - 2fj + fj-1)/h^2 + V fj  
%  = fj-1 (-1/2h^2) + fj (1/h^2 + V) + fj+1 (-1/2h^2)

clear; close all;

N=300;xmin=-10;xmax=10; h=(xmax-xmin)/N;   %cell centered with ghost
x=(xmin - h/2):h:(xmax + h/2);
x=x'; 

f = 0*x; %most of it is zero
wellleft = -1.5; wellright=1.5;
for j=1:N+2  %not the fastest way to do this, but it works
    if (x(j)<wellleft) | (x(j)>wellright) 
        f(j) = 5000;        
    end
end
        
c1 = -1/(2*h^2);
c2 = 1/h^2 + f;
c3 = -1/(2*h^2);
A = zeros(N+2,N+2);
B = eye(N+2,N+2);

for j=2:N+1         %rows 2 to N+1
    A(j,j-1) = c1;
    A(j,j)   = c2(j);
    A(j,j+1) = c3;
end

A(1,1)=.5; A(1,2)=.5;  B(1,1)=0; %first row
A(N+2,N+1)=.5; A(N+2,N+2)=.5; B(N+2,N+2)=0;  %last row

%%%%%start eigen.m code%%%%%%
[V,D]=eig(A,B);  % find the eigenvectors and eigenvalues

%w2raw=-g*diag(D);  % convert lambda to omega^2
w2raw = diag(D);

[epsilon,k]=sort(w2raw);  % sort omega^2 into ascending along with a
                     % sort key k(n) that remembers where each
                     % omega^2 came from so we can plot the proper
                     % eigenvector in V

for n=1:N  % run through the sorted list and plot each eigenvector
   % load the plot title into t
   gn=V(:,k(n)); % extract the eigenvector
   gnnormalized = gn/sqrt((h*sum(gn.*gn)));
   plot(x,gnnormalized,'b-');  % plot the eigenvector that goes with omega^2
   hold on
   plot(x,gnnormalized.^2,'r-'); % also plot psi^2
   plot(x,f/max(f),'g-');
   hold off
   t=sprintf(' epsilon = %g , integral of psi^2 = %g',epsilon(n),h*sum(gnnormalized.^2));
   title(t);xlabel('x');ylabel('g(n,x)');  % label the graph
   pause
end
%%%%%end eigen.m code%%%%%%



% first five eigenvalues from matlab: epsilon = .667882, 2.3929, 4.69426, 7.32979, 10.2329
% not certain they are correct... all my amplitudes were 1, not what manual has
% later: yes, they are correct... book is displaying *normalized* psi^2.

