%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% numerov.m
%
% solution of Schrodinger Equation for Cornell potential
% using numerov method. Program outputs eigen values only.
%
% July 2024
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;
tic
N    = 1000;
xmin = -0;
xmax = +20;
h = (xmax-xmin)/N;

K = 1; % m/hbar^2 value
tinyValue = 0.01;
toler = 1e-12;

E = 0;

for j = 1:15

dE = 0.1;

V = @(x)   -0.4./(x+eps) + 0.8*x;
%V = @(x)   x.*x/2;
F = @(x,E) 2*K*(E-V(x));
x = linspace(xmin,xmax,N+1);
c = h^2/12;

psi = zeros(1,N+1);
psi(2) = tinyValue;

for i=2:N
    psi(i+1) = (  2*(1-5*c*F(x(i),E))*psi(i) - (1+c*F(x(i-1),E))*psi(i-1) ) / ...
               (1 + c*F(x(i+1),E));
end
psiold = psi(end);


while 1
  psi = zeros(1,N+1);
  psi(2) = tinyValue;
  E = E + dE;
  for i=2:N
     psi(i+1) = (  2*(1-5*c*F(x(i),E))*psi(i) - (1+c*F(x(i-1),E))*psi(i-1) ) / ...
                (1 + c*F(x(i+1),E));
  end

%    plot(x,4*psi,'-b')
%    hold on
%    plot(x,V(x),'-r')
%    hold off
%    grid
%    axis([-10 10 -2 10])
%    pause(0.01)

  if psi(end)*psiold < 0
    E = E - dE;
    dE = dE / 2;

    if dE < toler
        fprintf('%f ',E)
        E = E + 0.1; % this for the next eigen value
        break
    end
    
  end

end

end
dt = toc
