[x,y] = meshgrid([-9.995:0.1:9.995],[-9.995:0.1:9.995]);
r = sqrt(x.^2 + y.^2);
% r = 0.01:0.01:20;
R20 = (1/sqrt(2))*(1 - 0.5*r).*exp(-0.5*r);
R21 = (1/sqrt(24))*r.*exp(-0.5*r);

psi0 = R20/sqrt(4*pi);
psi1 = R21*sqrt(3/(4*pi)).*x./r;
psi2 = R21*sqrt(3/(4*pi)).*y./r;

f0 = sqrt(2/3)*(psi1 + (1/sqrt(2))*psi0);
f1 = sqrt(2/3)*((sqrt(3)/2)*psi2 - 0.5*psi1 + (1/sqrt(2))*psi0);
f2 = sqrt(2/3)*(-(sqrt(3)/2)*psi2 - 0.5*psi1 + (1/sqrt(2))*psi0);

figure(1)
contour(f0.^2)
% axes([0 200],[0 200])
% axis('square')
% title('f0')
