<<"PhysicalConstants`"
<<"Units`"

n=1;
l=0;
m=0;
L=0.1;
phi=0;
rad[z_,x_]=(Sqrt[(x^2)+(z^2)]);
ang[z_,x_]=ArcTan[x/z];
En=Convert[-((ElectronCharge^2)/(4 Pi VacuumPermittivity))(1/( 2BohrRadius))((1/(n^2))),ElectronVolt]
β = Convert[Nano Meter(Sqrt[-((8 ElectronMass En)/(PlanckConstantReduced^2))]) ,1];
Psi[r_,theta_,phi_,n_,l_,m_]=Exp[-β r/2](((β r))^l)LaguerreL[n-l-1, 2l+1,β r]SphericalHarmonicY[l,m,theta,phi];
DensityPlot[((Abs[(Psi[rad[z,x],ang[z,x],phi,n,l,m])])^2),{x,-L,L},{z,-L,L},AspectRatio->1,Mesh->False,PlotPoints->70,PlotRange->All]

n=4;
l=0;
m=0;
L=1.5;
phi=0;
rad[z_,x_]=(Sqrt[(x^2)+(z^2)]);
ang[z_,x_]=ArcTan[x/z];
En=Convert[-((ElectronCharge^2)/(4 Pi VacuumPermittivity))(1/( 2BohrRadius))((1/(n^2))),ElectronVolt]
β = Convert[Nano Meter(Sqrt[-((8 ElectronMass En)/(PlanckConstantReduced^2))]) ,1];
Psi[r_,theta_,phi_,n_,l_,m_]=Exp[-β r/2](((β r))^l)LaguerreL[n-l-1, 2l+1,β r]SphericalHarmonicY[l,m,theta,phi];
Plot[((Abs[(Psi[rad[0,x],0,phi,n,l,m])])^2),{x,-L,L},AspectRatio->1,Mesh->False,PlotPoints->70,PlotRange->{0,0.1}]