(* Screened Coulomb potential -1/r + g (1-Exp[-delta r])/r *) (* Input: na, la - Principal and azimuthal quantum numbers g - Parameter of the potential nm - Order of perturbation theory Example of input line: na=2; la=1; g=.; nm=20; *) nq=na-1; (* Quantum number *) nm = morder+1; (* Incrementing morder because the last coefficient is calculated incorrectly because of undetected bug! *) nf=nm; kmi=0; kma=0; vf[1]=g; Do[vf[i]=-vf[i-1]/i,{i,2,nf}]; np=Floor[(nm+kmi)/2]+1; np=Min[np,na-la]; dn=na; ala=la*(la+1); mm=np+Floor[(nm+kma)/2]; dnh=dn/2; dnsq=dn^2/2; mm1=mm+1; Do[ Do[t[n,i]=0,{n,1,nm}]; a=dn+i-np; x1[i]=2*a; x2[i]=a*(a+1)-ala, {i,1,mm}]; energy[0]=-1/2/na^2; Print["E0 = ",energy[0]//InputForm," (Unperturbed Coulomb energy)"]; Do[ n9=n-1; ma1=Min[np+n,np+kma+nm-n]; mi1=Max[np-n,np-kmi-nm+n,1]; Do[ a=0; Do[x[i]=0,{i,1,mm1}]; x[nq]=1; If[n==1,Goto[22]]; Do[nt=n-nx; If[nx>nf,Goto[30]]; c=0; ma2=Min[np+nt,nq+nx]; mi2=Max[np-nt,nq-nx,1]; b=0; If[mi2>1,b=x[mi2-1]]; Do[a1=x2[nr]*x[nr+1]+x1[nr]*x[nr]+b; b=x[nr]; x[nr]=a1; c=c+a1*t[nt,nr], {nr,mi2,ma2}]; a=a-c*vf[nx]*dn*dnh^nx; Label[30];a1=x1[nq]*t[nt,nq]; If[nqmi2,a1=a1+x2[nq-1]*t[nt,nq-1]]; a=a+a1*e[nx]*dnsq, {nx,1,n9}]; Label[22];b=0; If[np>1,b=x[np-1]]; If[n",n," = ",energy[n]//InputForm]; If[npmi1,t[n,np-1]=t[n,np-1]-e[n]*dnsq]; t[n,np]=0, {n,1,nm-1(*!*)}];