(* Screened Coulomb potential -1/r + g (1-Exp[-delta r])/r *)
(* for symbolic g *)
(* the last coefficient is calculated incorrectly because of a bug! *)
nm=111;
Timing[
(*SetDirectory["~/math/screened"];*)
os=OpenWrite["series.dat"];
SetOptions[os,PageWidth->Infinity];
Clear[g];
nf=nm;
kmi=0;
kma=0;
na=1;
la=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}];
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;
   Clear[en];
   en[n]=e[n];
   Print[Definition[en]];
   Save[os,en];
   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}];
Close[os];
]//First