%------------------------------------------------------------
% Tony Hyun Kim
% 8.06 Term Paper on Hartree-Fock Approximation
%------------------------------------------------------------
% elec2_se: 
%   Numerically integrate the Schroedinger equation
% Parameters:
%   - r: the (radial) domain over which to solve SE
%   - V_ee: the potential
%------------------------------------------------------------
function psi = elec2_se(r,V)

% We are solving the radial equation that arises after 
% separating the 3D Hamiltonian. (See Griffiths, Ch.4. In
% particular, page 140)

% I let the angular momentum quantum number, l = 0, so that
% there is no centrifugal potential term.

% Note that we assume a spherically symmetric solution, so
% psi = R(r)*Y(theta,phi) = R(r). Since we define R = u(r)/r,
% where u(r) is the function appearing in the radial equation,
% the pertinent normalization condition is that integral from
% 0 to infty of |u(r)|^2 is 1/(4*pi) [Due to my convention, 
% this result differs from that given in Griffiths.]

% Furthermore, since R = u(r)/r, to keep R bounded, recall 
% that we must have u(r=0) = 0

% We will also fix d_psi(r=0) to be 1, since we will be 
% scaling it with a normalization coefficient 
% in the end, anyway
              
% Hence, the only variable available to us, is the energy that makes
% u bounded (i.e. lim r->infty u(r) = 0) without any internal nodes
% (the latter condition guarantees lowest energy state)

dr = r(2)-r(1);
N  = length(r);

% Convergence criterion
conv = 0.0001;
u_best = 1; % The variable where we store the results

% Initial guess
E       = linspace(-1,-0.6,100);

while(u_best > conv)
    u_end = zeros(1,length(E));
    for n = 1:length(E)
        % Reset psi and d_psi
        u   = zeros(1,N);
        d_u = ones(1,N);
    
        % In the inner loop, we integrate d^2u/dr^2 = -2(E-V)u
        for k = 2:N
            u(k)   = u(k-1)   + dr * d_u(k-1);
            d_u(k) = d_u(k-1) + dr * -2*(E(n)-V(k))*u(k);
        end
    
        u_end(n) = u(N);
    end

    % Evaluate results
    %plot(E,abs(u_end));
    %drawnow;
    [u_best best_index] = min(abs(u_end));
    E_best = E(best_index);
    % Next run
    E = linspace(E(best_index-1),E(best_index+1),100);
end

% Return the best result
u   = zeros(1,N);
d_u = ones(1,N);
for k = 2:N
    u(k)   = u(k-1)   + dr * d_u(k-1);
    d_u(k) = d_u(k-1) + dr * -2*(E_best-V(k))*u(k);
end
% Compute normalization, A
integ = 0;
for k = 1:N
    integ = integ + dr*u(k)^2;
end
A = 1/sqrt(4*pi*integ);
psi = A*u./r; 

