(* 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;
*)
(* Testing that na, la lie in a legitimate range *)
If[IntegerQ[la]&&la>=0&&IntegerQ[na]&&na>la,
Print["Principal quantum number: n = ",na];
Print["Azimuthal quantum number: l = ",la],
Print["Quantum numbers must be non-negative integers and n>l!"];
Exit[]];
(* Testing that g lies in a legitimate range *)
If[NumberQ[g],
   Print["Parameter g = ",g//InputForm],
   If[(!ValueQ[g]),
      Print["Parameter g has no value."],
      Print["Parameter g must be a number or has no value!"];
      Exit[]]];
(* Testing that nm lies in a legitimate range *)
If[IntegerQ[nm]&&nm>0,Print["Order of perturbation theory: N = ",nm],
Print["Order of perturbation theory must be a positive integer!"];
Exit[]];
nq=na-1; (* Quantum number *)
nm=nm+1(*!*);
(* 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}];
Print["Printing PT coefficients (expansion in delta) for screened Coulomb potential  -1/r + g (1-Exp[-delta r])/r ...\n"];
energy[0]=-1/2/na^2;
Print["E[0] = ",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[",n,"] = ",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(*!*)}];
(* The last coefficient is not calculated because of a bug in the program! *)
morder=nm;
If[(!ValueQ[g]),Exit[]];