function dydr=hydrogendy(r,y)
    
%  -1/2 f''(r) -1/r f(r) = Ef(r)
% (Hydrogen radial equation)
%        
    
    global E l ;

    V    = -1/r + l*(l+1)./r.^2 ;  % V + centrifugal
    dydr(1) = y(2)  ;        % this is f'(r)  = g(r) = y(2)
    dydr(2) = -2*(E-V)*y(1) ; % this is f''(r) = g'(r) = -2(E-V)*f(r) = -2(E-V)*y(1)
    dydr = dydr';