Psi[r_,theta_,phi_,n_,l_,m_]=Exp[- r/2]((( r))^l)LaguerreL[n-l-1, 2l+1, r]SphericalHarmonicY[l,m,theta,phi];
azimuth[x_,y_,z_]:=If[z==0,(π/2),ArcTan[((Sqrt[(x^2)+(y^2)])/z)]]
polar[x_,y_]:=If[x==0,(π/2),ArcTan[(y/x)]]
ψ[x_,y_,z_,n_,l_,m_]:=Simplify[((Abs[Psi[(Sqrt[(x^2)+(y^2)+(z^2)]),azimuth[x,y,z],polar[x,y],n,l,m]])^2),Assumptions->{x∈Reals,y∈Reals,z∈Reals}]

L=8;
ContourPlot[ψ[x,0,z,2,1,0],{x,-L,L},{z,-L,L}]

L=1;
ContourPlot3D[ψ[x,y,z,1,0,0]==1/8/π,{x, -L,L}, {y, -L, L}, {z, -L, L}]

L=3;
ContourPlot3D[{ψ[x,y,z,1,0,0]==1/8/π,ψ[x,y,z,1,0,0]==1/16/π,ψ[x,y,z,1,0,0]==1/32/π},{x, -L,L}, {y, -L, L}, {z, -L, L},Mesh->None,ContourStyle->{Directive[Blue,Opacity[0.2],Specularity[White,30]],Directive[LightGreen,Opacity[0.2],Specularity[White,30]],Directive[Orange,Opacity[0.2],Specularity[White,30]]}]

L=6.2;
ContourPlot3D[{ψ[x,y,z,2,1,0]==1/16/π,ψ[x,y,z,2,1,0]==1/8/π,ψ[x,y,z,2,1,0]==1/4/π},{x, -L,L}, {y, -L, L}, {z, -L, L},Mesh->None,ContourStyle->{Directive[Blue,Opacity[0.4],Specularity[White,30]],Directive[Blue,Opacity[0.3],Specularity[White,30]],Directive[Blue,Opacity[0.2],Specularity[White,30]]}]