(* 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["E<SUB>0</SUB> = ",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[nq<ma2,a1=a1+t[nt,nq+1]];
         If[nq>mi2,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<nf,a=a-(x2[np]*x[np+1]+x1[np]*x[np]+b)*vf[n]*
       dn*dnh^n];
      If[Not[nq==np],a=a/(nq-np)];
      t[n,nq]=a//Expand,
   {nq,mi1,ma1}];
   e[n]=(-t[n,np]/dn^3)//Expand;
   energy[n]=e[n];
   Print["E<SUB>",n,"</SUB> = ",energy[n]//InputForm];
   If[np<ma1,t[n,np+1]=t[n,np+1]+e[n]*x2[np]*dnsq];
   If[np>mi1,t[n,np-1]=t[n,np-1]-e[n]*dnsq];
   t[n,np]=0,
{n,1,nm-1(*!*)}];