% ---------------------------------
% hydrogen radial wave function
% ---------------------------------
clear all;close all;
N=300;
r0 = 0.1;  % smallest r
rinf = 30;  % largest r
H=zeros(N,N); % Hamiltonian matrix
r = linspace(r0,rinf,N);
h = r(2)-r(1); % stepsize  
alpha = 2.0; % Coulomb: alpha = 2.0 
for i=1:N
        H(i,i) = 2/h^2-alpha/r(i);
    for j=1:N 
        if abs(i-j)==1  
          H(i,j) = -1/h^2; 
        end
    end
end
[V,E]=eig(H);  % returns two matrices, E is diagonal
% sort eigenvalues & vectors, just to be sure
[E,ind]=sort(diag(E));
V = V(:,ind);
% output
for i=1:10
    if abs(alpha-2.0)<1e-10 
        fprintf('%5s %5d %20s %20.10f %20s %20.10f\n','i=',i,'E=',E(i),'exaxt E=',-1./i^2);
    else
        fprintf('%5s %5d %20s %20.10f\n','i=',i,'E=',E(i));
    end
end
plot(r,V(:,1),'b.');
