% --------------------------------- % 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.');