r=[1:1:1600];
ao=52.92;
Z=1;
A=(Z/ao)^(3/2);

p=2*Z*r/ao;
orb1s=A*2*exp(-p/2);

p=Z*r/ao;
orb2s=A*(1/(2*sqrt(2)))*(2-p).*exp(-p/2);

p=2*Z*r/(3*ao);
orb3s=A*(1/(9*sqrt(3)))*(6-6*p+p.^2 ).*exp(-p/2);

p=Z*r/ao;
orb2p=A*(1/(2*sqrt(6)))*p.*exp(-p/2);

p=2*Z*r/(3*ao);
orb3p=A*(1/(9*sqrt(6)))*(4-p).*exp(-p/2);

p=2*Z*r/(3*ao);
orb3d=A*(1/(9*sqrt(30)))*p.^2.*exp(-p/2);

term=4*pi*r.^2;
orb1s2=term.*orb1s.^2;
orb2s2=term.*orb2s.^2;
orb3s2=term.*orb3s.^2;
orb2p2=term.*orb2p.^2;
orb3p2=term.*orb3p.^2;
orb3d2=term.*orb3d.^2;

last=1000;
figure(1)
subplot(1,2,1)
plot(r(1:last),orb1s(1:last),'r',r(1:last),orb2s(1:last),'g')
title('1s and 2s')
xlabel('radius (angstroms) -->')
subplot(1,2,2)
plot(r(1:last),orb1s2(1:last),'r',r(1:last),orb2s2(1:last),'g')
title('1s and 2s - 4*pi*r2*psi-squared')
xlabel('radius (angstroms) -->')
subplot(1,2,2)

last=1500;
figure(2)
subplot(1,2,1)
plot(r(1:last),orb2s(1:last),'r',r(1:last),orb3s(1:last),'g')
title('2s and 3s')
xlabel('radius (angstroms) -->')
subplot(1,2,2)
plot(r(1:last),orb2s2(1:last),'r',r(1:last),orb3s2(1:last),'g')
title('2s and 3s - 4*pi*r2*psi-squared')
xlabel('radius (angstroms) -->')
subplot(1,2,2)

last=1200;
figure(3)
subplot(1,2,1)
plot(r(1:last),orb2p(1:last),'r',r(1:last),orb3p(1:last),'g')
title('2p and 3p')
xlabel('radius (angstroms) -->')
subplot(1,2,2)
plot(r(1:last),orb2p2(1:last),'r',r(1:last),orb3p2(1:last),'g')
title('2p and 3p - 4*pi*r2*psi-squared')
xlabel('radius (angstroms) -->')
subplot(1,2,2)

last=1600;
figure(4)
subplot(1,2,1)
plot(r(1:last),orb3p(1:last),'r',r(1:last),orb3d(1:last),'g')
title('3p and 3d')
xlabel('radius (angstroms) -->')
subplot(1,2,2)
plot(r(1:last),orb3p2(1:last),'r',r(1:last),orb3d2(1:last),'g')
title('3p and 3d - 4*pi*r2*psi-squared')
xlabel('radius (angstroms) -->')
subplot(1,2,2)