% ---------------------------------
% hydrogen radial wave function
% ---------------------------------
clear all;close all;
N=300;
H=zeros(N,N);
h= 0.2;
rr = 0.2:h:10000;
r = rr(1:N);
h = r(2)-r(1)
for i=1:N
        H(i,i) = 2/h^2-0.1/r(i);
    for j=i-1:2:i+1 
        if (j>0 && j<N+1) 
          H(i,j) =  -1/h^2; 
        end
    end
end
[V,E]=eig(H);
for i=1:10
    disp(['i= ' num2str(i) ' E= ' num2str(E(i,i)) ' exact E = ' num2str(-1./i^2)])
end
plot(r,V(:,1))