% *****************************************************
% Eigenvalue equation
% Solution from left and right to a meeting point
% change eigenvalue until logarithmic derivatives match
%
% Vesa Apaja May 2011
% *****************************************************

clear all;close all;
global E l ; % energy and angular momentum eigenvalues

l=0; 
acc  = 1e-15;
options = odeset('RelTol',1e-13,'AbsTol',1e-15);


% r-space to store solutions 
rinf = 40;  % wave function must be zero before rinf
rmeet = 2;  % meeting point             
step = 0.1; 
r = 1e-6:step:rinf; 
meet = rmeet/step;
n = size(r,2); 
r1 = r(1:meet);
ii = 0;
for i=n:-1:meet
    ii = ii + 1;
    r2(ii) = r(i);
end

% energy limits
E1 = -5;
E2 = 0;

diff1 = 1;
while(true)
    % try new eigenvalue, bisection method
    E = (E1+E2)/2 ;
    % solve DY from left to right and right to left
    % Obs: Psi(r,theta,phi)=R_l(r)/r Y_lm(theta,phi)
    % this solves function R_l(r)
    % initial values: R(0)=0    R'(0)=small 
    %                 R(rinf)=0 R'(rinf)=small  (sign irrelevant)
    % 
    [r1,f1]=ode113(@hydrogendy, r1, [0  1e-8], options);
    [r2,f2]=ode113(@hydrogendy, r2, [0 -1e-8], options);
    % logarithmic derivatives and their difference
    log1 = f1(end,2)/f1(end,1);
    log2 = f2(end,2)/f2(end,1);
    diff = log1-log2;
    
    fprintf(' E = %20.15f   diff = %20.15f \n',E,diff)
    % change E1 or E2 to bracket E for a smaller diff
    if(diff*diff1>0) 
        E1 = E;
        diff1 = diff;
    else
        E2 = E;
    end            
    if(E2-E1<acc) 
        break;        
    end
end
fprintf(' E = %20.15f diff = %20.15f \n',E,diff)
% scale solutions to match and max abs value to 1
f2(:,1) = f1(end,1)/f2(end,1)*f2(:,1)
scale = max(max(abs(f1(:,1))),max(abs(f2(:,1))));
f1(:,1) = f1(:,1)/scale;
f2(:,1) = f2(:,1)/scale;

plot(r1,f1(:,1),'r-'); 
hold on;
plot(r2,f2(:,1),'g-')

    
