a = 1; (*plots in units of Bohr radius*)

F[n_,l_,r_] := r^l  Exp[-r/2] (n + l)! LaguerreL[n-l-1, 2l +1, r];

Plot[F[1,0,r], {r,0,10}]

ψ[n_, l_, r_] := a^(-3/2) 2/(n^2) Sqrt[(n-1-l)!/( ( n+l)!)^3] F[n,l,2 r / (n a)];

Plot[ψ[1,0,r], {r,0,10}]

P[n_,l_,r_] := ψ[n,l,r] r ;

Plot[P[1,0,r], {r,0,10}]

NIntegrate[P[2,0,r]*P[2,0,r], {r,0,200}]

NIntegrate[P[2,0,r]*P[1,0,r], {r,0,50}]

NIntegrate[P[1,0,r]*r *P[2,1,r], {r,0,50}]

dipol1snp[n_] :=NIntegrate[P[1,0,r]*r *P[n,1,r], {r,0,50}]

dipol1snp[2]

Polariz[N_]:= Sum[ 2/3 *(dipol1snp[n])^2 / (1/2  -  1/(2 n^2)), {n,2,N}]

Table[ Polariz[n], {n,2,50}]

Pol = {2.9596211056347688`,3.3600727657907874`,3.49201416112381`,3.5525104453045393`,3.5855166251036534`,3.605583111050563`,3.618724273297262`,3.6278117043904627`,3.634363555506195`,3.6392460074683433`,3.642983477164115`,3.645908953920794`,3.648242315088023`,3.6501335850026995`,3.6516880189379717`,3.652981256483239`,3.654068802517756`,3.654992142123482`,3.655782790780401`,3.656465040740532`,3.657057862626028`,3.6575762469858555`,3.658032166857579`,3.6584352790453676`,3.658793442203938`,3.6591131044903187`,3.659399597033214`,3.6596573585096874`,3.6598901087231073`,3.660100984009285`,3.66029264377712`,3.66046735501184`,3.6606270598032196`,3.6607734296888426`,3.660907909675926`,3.661031754123813`,3.6611460561634117`,3.661251771951088`,3.6613497407685567`,3.661440701762711`,3.6615253079525565`,3.661604138001688`,3.6616777061547126`,3.661746470657849`,3.6618108409224037`,3.6618711836411926`,3.6619278280292806`,3.6619810703294773`,3.66203117769819`}

ListPlot[Pol]