
f[10,rho_]:=2;
f[20,r_]:=1/(2*Sqrt[2])*(2-r);
f[21,r_]:=1/(2*Sqrt[6])*r;
f[30,r_]:=1/(9*Sqrt[3])*(6-4*r+4/9*r^2);
f[31,r_]:=1/(9*Sqrt[6])*(4-2/3*r)*2/3*r;
f[32,r_]:=1/(9*Sqrt[30])*4/9*r^2;
f[40,r_]:=1/96*(24-18*r+3*r^2-1/8*r^3);
f[41,r_]:=1/(64*Sqrt[15])*r*(20-5*r+1/4*r^2);
f[42,r_]:=1/(384*Sqrt[5])*r^2*(6-1/2*r);
f[43,r_]:=1/(768*Sqrt[35])*r^3;
R[n_,l_,Z_,r_]:=f[n*10+l,Z*r]*Z^(3/2)*Exp[-Z*r/n];
			

P[n_,l_,Z_,r_]:=r^2*R[n,l,Z,r]^2;
			
<<PlotLegends`
Plot[{P[1,0,1,r],P[2,0,1,r],P[2,1,1,r]},{r,0,10},PlotRange->{Full,Full},AxesLabel->{"r","P"},
			PlotLegend->{Style["1s",16],Style["2s",16],Style["2p",16]},LegendPosition->{1.1,-0.4},LegendSize->{0.3,0.3},LegendShadow->None,
			PlotStyle->Thick,BaseStyle->{FontSize->16},PlotLabel->Style["Fonction de distribution radiale",16]]
Plot[{P[3,0,1,r],P[3,1,1,r],P[3,2,1,r]},{r,0,20},PlotRange->{Full,Full},AxesLabel->{"r","P"},
			PlotLegend->{Style["3s",16],Style["3p",16],Style["3d",16]},LegendPosition->{1.1,-0.4},LegendSize->{0.3,0.3},LegendShadow->None,
			PlotStyle->Thick,BaseStyle->{FontSize->16},PlotLabel->Style["Fonction de distribution radiale",16]]
Plot[{P[4,0,1,r],P[4,1,1,r],P[4,2,1,r],P[4,3,1,r]},{r,0,40},PlotRange->{Full,Full},AxesLabel->{"r","P"},
			PlotLegend->{Style["4s",16],Style["4p",16],Style["4d",16],Style["4f",16]},LegendPosition->{1.1,-0.4},LegendSize->{0.3,0.3},LegendShadow->None,
			PlotStyle->Thick,BaseStyle->{FontSize->16},PlotLabel->Style["Fonction de distribution radiale",16]]

Y[n_,l_,theta_,phi_]:=SphericalHarmonicY[n,l,theta,phi];
			
SphericalPlot3D[Abs[Y[0,0,theta,phi]]^2,{theta,0,Pi},{phi,0,2*Pi},PlotLabel->"1s"]
SphericalPlot3D[Abs[Y[1,0,theta,phi]]^2,{theta,0,Pi},{phi,0,2*Pi},PlotLabel->"2p"]
SphericalPlot3D[Abs[Y[2,0,theta,phi]]^2,{theta,0,Pi},{phi,0,2*Pi},PlotLabel->"3d"]
SphericalPlot3D[Abs[Y[2,1,theta,phi]]^2,{theta,0,Pi},{phi,0,2*Pi},PlotLabel->"3d"]
SphericalPlot3D[Abs[Y[2,2,theta,phi]]^2,{theta,0,Pi},{phi,0,2*Pi},PlotLabel->"3d"]
