#####################################################################
##PEKERIS: Save this file as  PEKERIS                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read PEKERIS<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
##with crucial help from Christoph Koutschan                          #
##zeilberg at math dot rutgers dot edu                               #
######################################################################
 
#started April 20, 2010: 
#that revisits 
#"Ground State of Two-Electron Atom"  Phys. Rev. 112 (1958) 1649-1658
 
print(`Created: April-May 2010`):
print(` This is PEKERIS `):
print(`It is a Maple package that accompanies the article `):
print(` The 1958 Pekeris-Accad-WEIZAC Ground-Breaking`):
print(` Collaboration that computed Ground States of Two-Electron Atoms (and its 2010 Redux) `):
print(``):
print(`by  Christoph Koutschan Doron Zeilberger.`):
print(`This article revisits Chaim Leib Pekeris' seminal paper`):
print(`"Ground State of Two-Electron Atom"  Phys. Rev. 112 (1958) 1649-1658`):
print(`and redoes the computations cleverly programmed `):
print(`(in machine language!), by Yigal Accad`):
print(`and performed in the pioneering WEIZAC computer.`):
print(``):
print(` That article is available from both authors' websites.`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
print(` A Mathematica package, written by Christoph Koutschan is also`):
print(`available from his website`):
print(``):

 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` ApplyDXi, ApplyDXiMono, ApplyMono , Chaim1,ChaimPC,ChaimPC3, `):
 print(` Comps,CompsA, ConOp, ConOpLin,`):
 print(` ConvertMonomial, ConvertOp, , COV, EV, FastDet `):
 print(` FindRoot, FindRoot1 `):
 print(` GenPol, GuessPol, Hafukh, `):
 print(` InterpolateCP, Kefel, Klmn, MakeEqs,MakeEqsOrtho, `):
 print(` MakeEqsPara, MatrixOrtho,MatrixPara, MatrixParaYigal, Pek2KC`):
 print(`  SimplifyMonomial1, Yafe, YafePlus, Zinn `):



else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(` Chaim, EnergyOrtho,EnergyPara,PaperOrtho, PaperPara `):
 print(` Pekeris1, Pekeris2,  `):


elif nops([args])=1 and op(1,[args])=ApplyDXi then
print(`ApplyDXi(F,R,X,D1,Oper,i): Given a linear partial differential `):
print(`operator Oper in the diff. operators D1[R[1]],..., D1[R[m]]`):
print(`with coefficients that are functions of X[1],..., X[m]`):
print(`left-multiplies it by differentiation with respect to X[i].`):
print(`For example, try:`):
print(`ApplyDXi([x^2+y^2,x*y],[u,v],[x,y],E,2*x*E[u]+y*E[v],1);`):

elif nops([args])=1 and op(1,[args])=ApplyDXiMono then
print(`ApplyDXiMono(F,R,X,D1,Oper,i): Given a linear partial differential `):
print(`operator Oper in the diff. operators D1[R[1]],..., D1[R[m]]`):
print(`with coefficients that are functions of X[1],..., X[m]`):
print(`that is a monomial in D1[R[1]],...,D1[R[m]]`):
print(`left-multiplies it by differentiation with respect to X[i]`):
print(`For example, try:`):
print(`ApplyDXiMono([x^2+y^2,x*y],[u,v],[x,y],E,2*x*E[u],1);`):

elif nops([args])=1 and op(1,[args])=ApplyMono then
print(`ApplyMono(F,R,X,D1,Mono): `):
print(`Given a transormation `):
print(`X[1]=F[1](R[1],..R[m]), X[2]=F[2](R[1],..R[m])`):
print(`expressing n X's in terms of m R's, where R is a list`):
print(`of independent variables and X is a list of dependent variables`):
print(`and D1 is a letter indicating that we denote diff. w.r.t. a variable`):
print(`v by D1[v], and Mono is a monomial differential operator`):
print(`in the X's and D1[X]'s converts it to a diff. operator`):
print(`in D1[R[1]], ..., D1[R1[m]] with coeff. that are a`):
print(`functions of X[1], ..., X[n]`):
print(`For example, try:`):
print(`ApplyMono([x^2+y^2,x*y],[u,v],[x,y],E,x*E[x]*E[y]);`):



elif nops([args])=1 and op(1,[args])=Chaim then
print(`Chaim(r1,r2,r12,n,D1,E,Z,MaxD)`):
print(` the two electron operator in n dimensions`):
print(`using the variables `):
print(`r1 (distance of electron ONE from the origin)`):
print(`r2 (distance of electron TWO from the origin)`):
print(`r12 (distance between the two electrons)`):
print(`D1[r1], D1[r2] , D1[r12] denote differentations`):
print(`w.r.t. the variables`):
print(`and E and Z are as in Pekeris paper (E is the energy and`):
print(`Z in the nuclear charge. For example, try:`):
print(`Chaim(r1,r2,r12,3,D1,E,Z,3);`):

elif nops([args])=1 and op(1,[args])=ChaimPC then
print(`ChaimPC(r1,r2,r12,n,D1,E,Z)`):
print(`The pre-computed Chaim operator for dimension n`):
print(`Do:`):
print(`ChaimPC(r1,r2,r12,3,D1,E,Z);`):

elif nops([args])=1 and op(1,[args])=ChaimPC3 then
print(`ChaimPC3(r1,r2,r12,D1,E,Z)`):
print(`The pre-computed Chaim operator for dimension 3`):
print(`Do:`):
print(`ChaimPC(r1,r2,r12,D1,E,Z);`):

elif nops([args])=1 and op(1,[args])=Chaim1 then
print(`Chaim1(x,y,n,D1,E,Z)`):
print(` the two electron operator in n dimensions`):
print(`using x[1], ..., x[n] for the coordinates of the first`):
print(`electron, y[1], ..., y[n] for the coordinates of the second`):
print(`electron and D1 to denote differentiation`):
print(`and E and Z are as in Pekeris' paper (E is the energy and`):
print(`Z in the nuclear charge. For example, try:`):
print(`Chaim1(x,y,3,D1,E,Z);`):


elif nops([args])=1 and op(1,[args])=Comps then
print(`Comps(d,k): all the vectors of length k that add-up to d exactly`):
print(`For example, try:`):
print(`Comps(3,2);`):

elif nops([args])=1 and op(1,[args])=CompsA then
print(`CompsA(d,k): all the vectors of length k that add-up to <=d`):
print(`For example, try:`):
print(`Comps(3,2);`):

elif nops([args])=1 and op(1,[args])=ConOp then
print(`ConOp(F,R,X,D1,Ope,MaxD): `):
print(`Given a transormation `):
print(`X[1]=F[1](R[1],..R[m]), X[2]=F[2](R[1],..R[m])`):
print(`expressing n X's in terms of m R's, where R is a list`):
print(`of independent variables and X is a list of dependent variables`):
print(`and D1 is a letter indicating that we denote diff. w.r.t. a variable`):
print(`v by D1[v], and Ope is a linear differential operator`):
print(`in the X's and D1[X]'s, converts Ope  to a diff. operator`):
print(`in D1[R[1]], ..., D1[R1[m]] and `):
print(`functions of R[1], ..., R[m]`):
print(`For example, try:`):
print(`ConOp([x^2+y^2,x*y],[u,v],[x,y],E,E[x]^2+E[y]^2,3);`):


elif nops([args])=1 and op(1,[args])=ConOpLin then
print(`ConOpLin(F,R,X,D1,Ope): `):
print(`Given an affine LINEAR transormation `):
print(`X[1]=F[1](R[1],..R[m]), X[2]=F[2](R[1],..R[m])`):
print(`expressing n X's in terms of m R's, where R is a list`):
print(`of independent variables and X is a list of dependent variables`):
print(`and D1 is a letter indicating that we denote diff. w.r.t. a variable`):
print(`v by D1[v], and Ope is a linear differential operator`):
print(`in the X's and D1[X]'s, converts Ope  to a diff. operator`):
print(`in D1[R[1]], ..., D1[R1[m]] and `):
print(`functions of R[1], ..., R[m]`):
print(`For example, try:`):
print(`ConOpLin([2*x+3*y+1,4*x+5*y-1],[u,v],[x,y],E,x*E[x]^2+y*E[y]^2);`):


elif nops([args])=1 and op(1,[args])=ConvertMonomial then
print(`ConvertMonomial(u,v,w,l,m,n,L,M,N,Li,Mono): Given`):
print(`a monomial in u,v,w, Mono, and discrete variables`):
print(`l,m,n, corresponding to them(resp.) and a symbol S1`):
print(`for the shift and a list L, of length 3, of operators representating`):
print(`left multiplication by u,v,w resp. finds the corresponding`):
print(`recurrence operator in terms of S1[l],S1[m],S1[n]`):
print(`For example try:`):
print(`ConvertMonomial(u,v,w,l,m,n,L,M,N,[-(l+1)*L+2*l+1-l*L^(-1),-(m+1)*M+2*m+1-m*M^(-1),-(n+1)*N+2*n+1-n*N^(-1)],u*v*w);`):


elif nops([args])=1 and op(1,[args])=ConvertOp then
print(`ConvertOp(u,v,w,l,m,n,L,M,N,Li,Ope): Given`):
print(`an operator in u,v,w, Ope, and discrete variables`):
print(`l,m,n, corresponding to them and respective shift`):
print(`operators L,M,N `):
print(`for the shift and a list Li, of length 3, of operators representating`):
print(`left multiplication by u,v,w resp. finds the corresponding`):
print(`recurrence operator in terms of L,M,N`):
print(`For example try:`):
print(`ConvertOp(u,v,w,l,m,n,L,M,N,[-(l+1)*L+2*l+1-l*L^(-1),-(m+1)*M+2*m+1-m*M^(-1),-(n+1)*N+2*n+1-n*N^(-1)],u+v+w);`):

elif nops([args])=1 and op(1,[args])=COV then
print(`COV(F,R,X,D1): `):
print(`Inputs a list F of functions (of size m, say)`):
print(`in the list of variables X (of length n, say)`):
print(`and a list of variables R (of length m)`):
print(`and a symbol D1 for differentiation`):
print(`outputs the list of first-order differential operators`):
print(`[D[X[1]],D[X[2]], ... ,D[X[n]]]`):
print(`as linear combinations of the D[R[1]], ..., D[R[m]]`):
print(`using the multivariable chain rule`):
print(`For example, try:`):
print(`COV([x^2+y^2,x*y],[u,v],[x,y],E);`):
print(``):


elif nops([args])=1 and op(1,[args])=EnergyOrtho then
print(`EnergyOrtho(omega,Z,d,memad,orekh): In memad dimensions`):
print(`inputs a value for omega >=5`):
print(`and positive integers d  and orekh such that`):
print(`omega-(orekh-1)*d>=3 and orekh>=5 `):
print(`(the truncation parameter in Pekeris' pertrubation method`):
print(`described in his above-mentioned article)`):
print(`outputs, under the ortho (i.e. anti-symmetric) assumption,`):
print(` estimates for N-thlargest eignevalue (E=-e^2)`):
print(`in atomic units, `):
print(`for the two-electron charge with charge Z`):
print(`using omega-(orekh-1)*d,omega-(orekh-2)*d,... , omega-d,omega`):
print(`and using interpolation`):
print(`The output is the best estimate  `):
print(`followed by the interpolated estimate`):
print(`for the largest eigenvalue `):
print(`For example, for the  for the square-root of the `):
print(`negative of the ground state energy`):
print(`for the 2 {}^3S Stae of Helium given in Table V`):
print(`of Pekeris' 1959 paper (Phys. Rev. v. 115, Nu. 5, 1216-1221`):
print(` page 1220 `):
print(`type :`):
print(`EnergyOrtho(19,2,3,3,5);`):

elif nops([args])=1 and op(1,[args])=EnergyPara then
print(`EnergyPara(omega,Z,d,memad,orekh): In memad dimensions`):
print(`inputs a value for omega >=5`):
print(`and positive integers d  and orekh such that`):
print(`omega-(orekh-1)*d>=3 and orekh>=5 `):
print(`(the truncation parameter in Pekeris' pertrubation method`):
print(`described in his above-mentioned article)`):
print(`outputs, under the para (i.e. symmetric) assumption,`):
print(` estimates for highest eigenvalue e (E=-e^2)`):
print(`in atomic units, `):
print(`for the two-electron charge with charge Z`):
print(`using omega-(orekh-1)*d,omega-(orekh-2)*d,... , omega-d,omega`):
print(`and using interpolation`):
print(`The output is the best estimate  `):
print(`followed by the interpolated estimate`):
print(`for the N-th largest eigenvalue `):
print(`For example, for the  for the square-root of the`):
print(`negative of the ground state energy`):
print(`of Helium, given in the Z=2 line in Table III of Pekeris' 1958 paper`):
print(`type :`):
print(`EnergyPara(11,2,1,3,5);`):

elif nops([args])=1 and op(1,[args])=EV then
print(`EV(C,e): Given a square-matrix C whose entries are affine-linear`):
print(`expressions in e, finds the first N "generalized" eigenvalues.`):
print(`in floating-point approximations.`):
print(`For example, try:`):
print(`EV(MatrixPara(2,e,5,3),e);`):

elif nops([args])=1 and op(1,[args])=FastDet then
print(`FastDet(lu): computes a symbolic determint of`):
print(` a matrix lu with entries that depend on a single variable `):
print(`For example, try:`):
print(`FastDet([[e,1],[1,e]]);`):

elif nops([args])=1 and op(1,[args])=FindRoot then
print(`FindRoot(A,e,L,prec):`):
print(`inputs a matrix A whose entries are expressions`):
print(`in e and that has an eigenvalue in the interval L (of the form`):
print(`[a,b]), and the required precision`):
print(`outputs the interval [a,(a+b)/2] or [(a+b),b] that`):
print(`has it. Try:`):
print(`FindRoot(MatrixPara(2,e,5,3),e,[17/10,18/10],10);`):

elif nops([args])=1 and op(1,[args])=FindRoot1 then
print(`FindRoot1(A,e,L): inputs a matrix A whose entries are expressions`):
print(`and e and that has an eigenvalue in the interval L (of the form`):
print(`[a,b]), outputs the interval [a,(a+b)/2] or [(a+b),b] that`):
print(`has it. Try:`):
print(`FindRoot1(MatrixPara(2,e,5,3),e,[17/10,18/10]);`):

elif nops([args])=1 and op(1,[args])=GenPol then
print(`GenPol(X,d,a): a generic polynomial in the list of variables X`):
print(`of degree d with coefficients indexed by the symbol a,`):
print(`followed by the set of coefficients`):
print(`For example, try:`):
print(`GenPol([x,y],3,a);`):

elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(f,X,F,R,MaxD): Given an expression f in the variables`):
print(`X[1], ..., X[n] that is known to be a polynomial in`):
print(`R[1],...,R[m] with R[i] substituted for F[i] where`):
print(`the F[i] is an expression in the X[1],...,X[n]`):
print(`tries to find that polynomial in R[1],..., R[m]`):
print(`For example, try:`):
print(`GuessPol(sqrt(x^2+y^2)^3,[x,y],[sqrt(x^2+y^2)],[r],5);`):

elif nops([args])=1 and op(1,[args])=Hafukh then
print(`Hafukh(F,X,R): Given an affine linear transformation`):
print(`from a list of variables X to a list of variables R`):
print(`given by F, returns the inverse transformation.`):
print(`For example, try:`):
print(`Hafukh([x+3*y,2*x-y],[x,y],[u,v]);`):

elif nops([args])=1 and op(1,[args])=InterpolateCP then
print(`InterpolateCP(L): given a list L of 4 consecutive values `):
print(`guesses the limiting value using eqs. (33) and (34) in`):
print(`Pekeris' paper`):

elif nops([args])=1 and op(1,[args])=Kefel then
print(`Kefel(ope1,ope2,n,N): multiplies ope1(n,N) by ope2(n,N)`):
print(`For example, try:`):
print(`Kefel(N,n,n,N);`):

elif nops([args])=1 and op(1,[args])=Klmn then
print(`The linearization of the vectors (l,m,n) with l<=m. `):
print(`For example, try: Klmn(0,1,1);`):

elif nops([args])=1 and op(1,[args])=MakeEqs then
print(`MakeEqs(A,Z,e,om,memad): makes the set of equations for`):
print(`Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om`):
print(`in memad dimensions`):
print(`The output is the set of equations followed by the`):
print(`set of variables.For example, do:`):
print(`MakeEqs(A,Z,e,1,3);`):

elif nops([args])=1 and op(1,[args])=MakeEqsOrtho then
print(`MakeEqsOrtho(A,Z,e,om,memad): In mamad dimensions`):
print(`makes the set of equations for`):
print(`Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om`):
print(`under the ortho assumption (Eq. 25 in Pekeris' paper)`):
print(`i.e. that A(l,m,n)=-A(m,l,n) .`):
print(`The output is the set of equations followed by the`):
print(`set of variables. For example, do:`):
print(`MakeEqsOrtho(A,Z,e,1,3);`):

elif nops([args])=1 and op(1,[args])=MakeEqsPara then
print(`MakeEqsPara(A,Z,e,om,memad): makes the set of equations for`):
print(`Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om`):
print(`in memad dimensions, `):
print(`under the para assumption (Eq. 25 in Pekeris' paper)`):
print(`i.e. that A(l,m,n)=A(m,l,n) .`):
print(`The output is the set of equations followed by the`):
print(`set of variables.For example, do:`):
print(`MakeEqsPara(A,Z,e,1,3);`):

elif nops([args])=1 and op(1,[args])=MatrixOrtho then
print(`MatrixOrtho(Z,e,om,memad): In memad dimensions`):
print(`makes the matrix of`):
print(`Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om,`):
print(`with the ortho assumption (i.e. A(l,m,n)=-A(m,l,n))`):
print(`The output is the set of equations followed by the`):
print(`set of variables.For example, do:`):
print(`MatrixOrtho(Z,e,1,3);`):

elif nops([args])=1 and op(1,[args])=MatrixPara then
print(`MatrixPara(Z,e,om,memad): In memad dimensions`):
print(`makes the matrix of`):
print(`Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om,`):
print(`with the para assumption (i.e. A(l,m,n)=A(m,l,n))`):
print(`The output is the matrix expressed as a list of lists`):
print(`For example, try:`):
print(`MatrixPara(Z,e,1,3);`):

elif nops([args])=1 and op(1,[args])=MatrixParaYigal then
print(`MatrixParaYigal(Z,e,om,memad): In memad dimensions`):
print(`makes the matrix of`):
print(`Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om,`):
print(`exactly the way it was done in Pekeris' table`):
print(`with the para assumption (i.e. A(l,m,n)=A(m,l,n))`):
print(`The output is the matrix expressed as a list of lists`):
print(`For example, to get Table II of Pekeris' paper, type:`):
print(`MatrixParaYigal(Z,e,2,3);`):


elif nops([args])=1 and op(1,[args])=PaperOrtho then
print(`PaperOrtho(psi,x,y,r,u,v,w,F,MaxZ,omega,A,E,e,Z,l,m,n,d,orekh):`):
print(`In memad dimensions, outputs a paper about the ground energy`):
print(`psi,x,y,r,u,v,w,F,A,E,e,Z,m,l,n are symbols as in Pekeris'`):
print(`paper, MaxZ is the the maximal charge (i.e. Z is taken from 1 to MaxZ)`):
print(`For example, to reproduce`):
print(`Pekeris' 1959 paper (and more) type:`):
print(`PaperOrtho(psi,x,y,r,u,v,w,F,10,22,A,E,e,Z,l,m,n,3,5);`)

elif nops([args])=1 and op(1,[args])=PaperPara then
print(`PaperPara(psi,x,y,r,u,v,w,F,MaxZ,omega,A,E,e,Z,l,m,n,d,orekh):`):
print(`In memad dimensions, outputs a paper about the ground energy`):
print(`psi,x,y,r,u,v,w,F,A,E,e,Z,m,l,n are symbols as in Pekeris'`):
print(`paper, MaxZ is the the maximal charge (i.e. Z is taken from 1 to MaxZ)`):
print(`For example, to reproduce`):
print(`Pekeris' paper type:`):
print(`PaperPara(psi,x,y,r,u,v,w,F,10,11,A,E,e,Z,l,m,n,1,5);`)



elif nops([args])=1 and op(1,[args])=Pekeris1 then
print(`Pekeris1(u,v,w,D1,n): The diff. operator Chaim in terms of the `):
print(`perimetric coordinates u,v,w, in n dimensions`):
print(`and D1 denotes differentiation. Do:`):
print(`Pekeris1(u,v,w,D1,e,Z,3);`);

elif nops([args])=1 and op(1,[args])=Pekeris2 then
print(`Pekeris2(l,m,n,e,Z,L,M,N,n): The differnce`):
print(`operator corresponding to Pekeris1, in n dimensions`):
print(`and using shift operators L,M,N, corresponding`):
print(`to l,m,n `):
print(`Do:`):
print(`Pekeris2(l,m,n,e,Z,L,M,N,3);`):


elif nops([args])=1 and op(1,[args])=Pekeris2CK then
print(`Pekeris2CK(l,m,n,e,Z,L,M,N): pre-computed`):
print(`Pekeris2 using Christoph Koutschan's package`):
print(`Pekeris2CK(l,m,n,e,Z,L,M,N);`):

elif nops([args])=1 and op(1,[args])=SimplifyMonomial1 then
print(`SimplifyMonomial1(L,x,D,n,N): Given a list L of`):
print(`Rewriting Rules, where L[i] is an expression`):
print(`for x^(i-1)*D^*(i-1) as an operator in n, N, `):
print(`(up to nops(L))`):
print(`finds the simplified form an operator in n,N`):
print(`only for the Monomial of the x^a*D^b with a>=b>=nops(L)-1.`):
print(`For example, try:`):
print(`SimplifyMonomial1([-(n+1)*N+(2*n+1)-n/N,n-n/N,(n^2-n)*(N-2/N+1/N^2)],x,D,n,N,x^3*D^2);`):

elif nops([args])=1 and op(1,[args])=Yafe then
print(`Yafe(P,X): given a multivariable polynomial P`):
print(`in X[1],..,X[n] (X is a list of variables)`):
print(`collects and simplifies.`):
print(`For example, try: `):
print(`Yafe(x^2*z+2*x*y*z+y^2*z,[z]); `):

elif nops([args])=1 and op(1,[args])=YafePlus then
print(`YafePlus(P,Y,X,F,R): given a multivariable polynomial P`):
print(`in Y[1],..,Y[k] (Y is a list of variables)`):
print(`collects and simplifies`):
print(`where the coeffs. are functions of X[1], ..., X[n]`):
print(`believed to be expressible as rational functions`):
print(`in the R[1], .., R[k] where R[1]=F[1], ...`):
print(`for example, try`):
print(`YafePlus(x^2*z+2*x*y*z+y^2*z,[z],[x,y],[x+y],[r],2);`):


elif nops([args])=1 and op(1,[args])=Zinn then
print(`Zinn(resh,n): The method of Zinn-Justin. Inputs a list resh`):
print(`believed to be appx. of the form mu^n*n^theta`):
print(`returns the estimated theta followed by the mu.`):
print(`For example, try:`):
print(`Zinn([seq(1.1^i*(i+1),i=1..20)],19);`):

else
print(`There is no ezra for`,args):
fi:
 
end:


####CKs computation
Pekeris2CK:=proc(l,m,n,e,Z,L,M,N) local ck:
ck:=
-10-4*M*m*L+12*N*e*l-4*l/L^2*Z+24*e*M*m^2-4*M*L*l+4*l^2/L^2*Z+4*N*n*m+4*n*l/L+24*m^2/M*e-4*e*M^2*m^2+4*Z*L^2*l^2-16*m*Z*L-4*m^2/M^2*e+30*n*e*L-32*e+8*n^2/N*e-8*n*e*
L^2+8*M*Z*L+8*m^2*e*L+8*m*L*l+8*M*m*l-2*N*M*m+6*L-14*n+32*Z+4*n*m/M+4*n*L*l-16*Z*L*l^2+32*m*Z*l+28*e*L-40*m^2*e-2/N/L*n*e*l+36*m*e*L+6*M+8*N-12*l-12*m-16*l^2/L*Z+4*
n*M*m+4*N*n*l-2*n/N*M-4*M*l/L+8*N*e*n^2+4*n/N*m-12*e*M^2*m-8*M^2*e*l+4*m/M^2*e-2*n/N*L+16*N*n*e-8*M*e*L-2*N*n*M+52*e*L*l+8*m*l/L-2*N*m/M+8*N*m^2/M*n*e+4*n/N*l-4*m/M
*L+4*n/N*Z*M*m-72*n*e*l+4*l*e/L-4*N*l^2*e+8*m/M*l-4*N*e*L-96*m*e*l-2*N*n*L+24*l^2*Z-16*M*Z*l+36*M*e*l+4*N*n*m*Z/M+8*Z*M^2+12*Z*L^2*l-4*m^2/M*l*e/L+2*N*n*l*e/L-4*M*m
^2*l*e/L+4*m/M^2*e*l+6*l/L+6*m/M+4*n/N-64*e*l+4*n/N*e*l+6*n/N*e*L+6*n*l*e/L+54*n*e*L*l+20*N*n*e*l-6*N*n*e*L-4*N*e*L*l-2*n/N*m/M-2*n/N*M*m-2*N*n*m/M-2*N*n*M*m-48*n*e
+8*N*e-16*n^2*e-6*n^2+6*M*m+14*n/N*e*L*l-6*N*n*e*L*l-8*n*m+4*n*M+12*N*n+4*N*m-2*N*M-n^2/N^2+n/N^2+4*n^2/N+4*N*n^2-4*m/M*l/L-4*m/M*L*l-4*M*m*l/L-4*M*m*L*l-2*n/N*l/L-\
2*n/N*L*l-2*N*n*l/L-2*N*n*L*l-16*m/M*Z*l+8*m/M*Z*L-16*m*l*Z/L-16*m*Z*L*l-16*M*m*Z*l+8*M*m*Z*L+8*M*l*Z/L+8*M*Z*L*l-8*n/N*Z*l+4*n/N*Z*L-8*n*l*Z/L-8*n*Z*L*l-8*N*n*Z*l+
4*N*n*Z*L+4*N*l*Z/L+4*N*Z*L*l-8*n/N*Z*m+4*n/N*Z*M-8*n*m*Z/M-8*n*Z*M*m-8*N*n*Z*m+4*N*n*Z*M+4*N*m*Z/M+4*N*Z*M*m+4*n/N*e*m-N^2*n^2-3*N^2*n-2*N^2+6*L*l-16*m*l+8*m*L+8*M
*l-4*M*L-8*n*l+4*n*L+4*N*l-2*N*L+48*Z*m-2*N*l/L-2*N*L*l-12*m*Z/M-28*Z*M-44*Z*M*m+4*m^2/M^2*Z-4*m/M^2*Z-16*m^2/M*Z+24*m^2*Z-16*Z*M*m^2+4*Z*M^2*m^2+12*Z*M^2*m-12*l*Z/
L-28*Z*L-44*Z*L*l+48*Z*l+8*Z*L^2+16*n*Z-8*N*Z+8*m/M*l*Z/L+8*m/M*Z*L*l+8*M*m*l*Z/L+8*M*m*Z*L*l-8*n/N*Z+16*n*Z*l-8*n*Z*L-8*N*n*Z-8*N*Z*l+4*N*Z*L+16*n*Z*m-8*n*Z*M-8*N*
Z*m+4*N*Z*M+4*m*e/M+28*e*M+52*e*M*m-72*n*e*m+30*n*e*M+12*N*e*m-4*N*e*M+4*n/N*l*Z/L+4*n/N*Z*L*l+4*N*n*l*Z/L+4*N*n*Z*L*l-64*e*m+4*n/N*m*Z/M+4*N*n*Z*M*m+14*n/N*e*M*m+6
*n/N*e*M+6*n*m*e/M+54*n*e*M*m+20*N*n*e*m-6*N*n*e*M-4*N*e*M*m-2*n^2/N^2*e*L+2*n/N^2*e*L+8*n^2/N*e*l+2*n^2/N*e*L-16*n^2*e*l+2*n^2*e*L+12*m/M*e*l-4*m/M*e*L+12*m*l*e/L+
60*m*e*L*l-8*e*M^2+2*N*n*m*e/M-6*N*n*e*M*m-8*e*L^2-4*l^2/L^2*e+4*l/L^2*e+24*l^2/L*e-40*l^2*e+24*e*L*l^2-4*e*L^2*l^2-12*e*L^2*l-2*n^2/N^2*e*L*l+2*n/N^2*e*L*l-2*n^2/N
*l*e/L+2*n^2/N*e*L*l+2*n^2*l*e/L+2*n^2*e*L*l-8*m/M*e*L*l+60*M*m*e*l-12*M*m*e*L-4*M*l*e/L-12*M*e*L*l+8*n^2/N*e*m+2*n^2/N*e*M+2*n^2*m*e/M+2*n^2*e*M*m+8*N*n^2*e*m-2*N*
n^2*e*M-4*m^2/M^2*n*e+4*m/M^2*n*e+24*m^2/M*n*e-40*m^2*n*e+24*n*e*M*m^2-8*M*m*l*e/L-16*M*m*e*L*l-2/N/M*n*e*m-16*n^2*e*m+2*n^2*e*M-8*n*e*M^2-2*n^2/N*m*e/M+2*n^2/N*e*M
*m+2*N*n^2*m*e/M-2*N*n^2*e*M*m-4*n*e*M^2*m^2-12*n*e*M^2*m-4*l^2/L^2*n*e+4*l/L^2*n*e+24*l^2/L*n*e-40*l^2*n*e+24*n*e*L*l^2-4*n*e*L^2*l^2-12*n*e*L^2*l+8*N*n^2*e*l-2*N*
n^2*e*L+2*N*n^2*l*e/L-8*m*e*L^2+8*M*l^2*e-2*N*n^2*e*L*l-4*l^2/L^2*m*e+4*l/L^2*m*e+24*l^2/L*m*e-40*l^2*m*e+24*m*e*L*l^2-4*m*e*L^2*l^2-12*m*e*L^2*l+8*m/M*l^2*e+8*M*l^
2*m*e-4*M*l^2/L*e-4*M*e*L*l^2-4*m^2/M^2*e*l-4*m/M*l^2/L*e-4*m/M*e*L*l^2-4*M*l^2/L*m*e-4*M*m*e*L*l^2+24*m^2/M*e*l-40*m^2*e*l+24*M*m^2*e*l-4*M^2*m^2*e*l-12*M^2*m*e*l-\
4*m^2/M*e*L+8*m^2*l*e/L+8*m^2*e*L*l-4*M*m^2*e*L-32*n*m*e*l+8*n*m*e*L+16*N*m*e*l-4*N*m*e*L+8*n*M*e*l-4*m^2/M*e*L*l-4*M*m^2*e*L*l-4*n/N*m*l*e/L+16*n/N*m*e*l-4*n/N*m*e
*L*l-4*n/N*m*e*L+8*n*m*l*e/L+8*n*m*e*L*l-4*N*n*m*l*e/L+16*N*n*m*e*l-4*N*n*m*e*L*l-4*N*n*m*e*L-4*N*m*l*e/L-4*N*m*e*L*l-4*n/N*m/M*e*l-4*n/N*M*m*e*l-4*n/N*M*e*l+8*n*m/
M*e*l+8*n*M*m*e*l-4*N*M*e*l-2*n^2/N^2*e*M+2*n/N^2*e*M-4*N*n*m/M*e*l-4*N*n*M*m*e*l-4*N*n*M*e*l-4*N*m/M*e*l-4*N*M*m*e*l-4*N*m^2*e-2*n^2/N^2*e*M*m+2*n/N^2*e*M*m-8*n/N*
e*M^2-4*N*m^2*n*e-4*N*m^2/M^2*e+4*N*m/M^2*e+8*N*m^2/M*e-4*N^2*m*e/M-4*m^2*n/N*e+8*n/N*e*M*m^2-4*n/N*e*M^2*m^2-12*n/N*e*M^2*m-4*N*m^2/M^2*n*e+4*N*m/M^2*n*e-2*N^2*n^2
*m*e/M-6*N^2*n*m*e/M-4*l^2*n/N*e+8*n/N*e*L*l^2-8*n/N*e*L^2-4*N*l^2*n*e-4*N*l^2/L^2*e+4*N*l/L^2*e+8*N*l^2/L*e-4*N^2*l*e/L-4*n/N*e*L^2*l^2-12*n/N*e*L^2*l-4*N*l^2/L^2*
n*e+4*N*l/L^2*n*e+8*N*l^2/L*n*e-2*N^2*n^2*l*e/L-6*N^2*n*l*e/L
:
ck:
Yafe(ck,[L,M,N]):
end:
####End CKs computation

#COV(F,R,X,D1):
#Inputs a list F of functions (of size m, say)
#in the list of variables X (of length n, say)
#and a list of variables R (of length m)
#and a symbol D1 for differentiation
#outputs the list of first-order differential operators
#[D[X[1]],D[X[2]], ... ,D[X[n]]]
#as linear combinations of the D[R[1]], ..., D[R[m]]
#with coefficients that are functions of X
#using the multivariable chain rule
#For example, try:
#COV([x^2+y^2,x*y],[u,v],[x,y],E);
COV:=proc(F,R,X,D1) local i,j:
option remember:
if nops(F)<>nops(R) then
RETURN(FAIL):
fi:
[seq(add(normal(diff(F[j],X[i]))*D1[R[j]],j=1..nops(F)),i=1..nops(X))]:
end:





#ApplyDXiMono(F,R,X,D1,Oper,i): Given a linear partial differential 
#operator Oper in the diff. operators D1[R[1]],..., D1[R[m]]
#with coefficients that are functions of X[1],..., X[m]
#that is a monomial in D1[R[1]],...,D1[R[m]]
#left-multiplies it by differentiation with respect to X[i]
#For example, try:
#ApplyDXiMono([x^2+y^2,x*y],[u,v],[x,y],E,2*x*E[u],1);
ApplyDXiMono:=proc(F,R,X,D1,Oper,i) local gu,Deg1,i1,Oper1:
gu:=COV(F,R,X,D1):
Deg1:=[seq(degree(Oper,D1[R[i1]]),i1=1..nops(R))]:
Oper1:=normal(Oper/mul(D1[R[i1]]^Deg1[i1],i1=1..nops(R))):

if {seq(normal(diff(Oper1,D1[R[i1]])),i1=1..nops(R))}<>{0} then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:
expand(diff(Oper1,X[i])*mul(D1[R[i1]]^Deg1[i1],i1=1..nops(R))+
gu[i]*Oper):

end:




#ApplyDXi(F,R,X,D1,Oper,i): Given a linear partial differential 
#operator Oper in the diff. operators D1[R[1]],..., D1[R[m]]
#with coefficients that are functions of X[1],..., X[m]
#left-multiplies it by differentiation with respect to X[i]
#For example, try:
#ApplyDXi([x^2+y^2,x*y],[u,v],[x,y],E,2*x*E[u]+y*E[v],1);
ApplyDXi:=proc(F,R,X,D1,Oper,i) local j:

if not type(Oper,`+`) then
 RETURN( ApplyDXiMono(F,R,X,D1,Oper,i)):
fi:

add(ApplyDXiMono(F,R,X,D1,op(j,Oper),i),j=1..nops(Oper)):
end:


#ApplyMono(F,R,X,D1,Mono): 
#Given a transormation 
#X[1]=F[1](R[1],..R[m]), X[2]=F[2](R[1],..R[m])
#expressing n X's in terms of m R's, where R is a list
#of independent variables and X is a list of dependent variables
#and D1 is a letter indicating that we denote diff. w.r.t. a variable
#v by D1[v], and Mono is a monomial differential operator
#in the X's and D1[X]'s converts it to a diff. operator
#in D1[R[1]], ..., D1[R1[m]] with coeff. that are 
#functions of X[1], ..., X[n]
#For example, try:
#ApplyMono([x^2+y^2,x*y],[u,v],[x,y],E, x*E[x]*E[y]);
ApplyMono:=proc(F,R,X,D1,Mono) local i,j,L,lu,gu:

L:=[seq(degree(Mono,D1[X[i]]),i=1..nops(X))]:

lu:=normal(Mono/mul(D1[X[i]]^L[i],i=1..nops(X))):

if {seq(expand(normal(diff(lu,D1[X[i]]))),i=1..nops(X))}<>{0} then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

gu:=1:

for i from 1 to nops(L) do
for j from 1 to L[i] do:
gu:=ApplyDXi(F,R,X,D1,gu,i):
od:
od:
expand(lu*gu):
end:








#Yafe(P,X): given a multivariable polynomial P
#in X[1],..,X[n] (X is a list of variables)
#collects and simplifies
#for example, try
#Yafe(x^2*z+2*x*y*z+y^2*z,[z]);
Yafe:=proc(P,X) local i,lu,P1,v,gu,i1,gu1,L:

if not type(P,`+`) then
 RETURN(normal(P)):
fi:

lu:={}:

for i from 1 to nops(P) do
P1:=op(i,P):
v:=[seq(degree(P1,X[i1]),i1=1..nops(X))]:
lu:=lu union {v}:
od:

gu:=0:

for L in lu do
gu1:=P:

for i1 from 1 to nops(X) do
 gu1:=coeff(gu1,X[i1],L[i1]):
od:

gu:=gu+ factor(normal(gu1))*mul(X[i1]^L[i1],i1=1..nops(X)):
od:

gu:

end:




#Chaim1(x,y,n,D1,E,Z)
# the two electron operator in n dimensions
#using x[1], ..., x[n] for the coordinates of the first
#electron, y[1], ..., y[n] for the coordinates of the second
#electron and D1 to denote differentiation
#and E and Z are as in Pekeris' paper (E is the energy and
#Z in the nuclear charge. For example, try:
#Chaim1(x,y,3,D1,E,Z);
Chaim1:=proc(x,y,n,D1,E,Z) local r1,r2,r12,i:
r1:=sqrt(add(x[i]^2,i=1..n)):
r2:=sqrt(add(y[i]^2,i=1..n)):
r12:=sqrt(add((x[i]-y[i])^2,i=1..n)):
add(D1[x[i]]^2,i=1..n)+add(D1[y[i]]^2,i=1..n)+
2*(E+Z/r1+Z/r2-1/r12):
end:


#Chaim1a(x,y,n,D1,E,Z,r)
# the two electron operator in n dimensions
#using x[1], ..., x[n] for the coordinates of the first
#electron, y[1], ..., y[n] for the coordinates of the second
#electron and D1 to denote differentiation
#and E and Z are as in Pekeris' paper (E is the energy and
#Z in the nuclear charge. For example, try:
#Chaim1a(x,y,3,D1,E,Z,r);
Chaim1a:=proc(x,y,n,D1,E,Z,r) local r1,r2,r12,i:
add(D1[x[i]]^2,i=1..n)+add(D1[y[i]]^2,i=1..n)+
2*(E+Z/r1+Z/r2-1/r12):
end:





#Chaim(r1,r2,r12,n,D1,E,Z,MaxD)
# the two electron operator in n dimensions
#using the variables 
#r1 (distance of electron ONE from the origin)
#r2 (distance of electron TWO from the origin)
#r12 (distance between the two electrons)
#D1[r1], D1[r2] , D1[r12] denote differentations
#w.r.t. the variables
#and E and Z are as in Pekeris' paper (E is the energy and
#Z in the nuclear charge. For example, try:
#Chaim(r1,r2,r12,3,D1,E,Z,3);
Chaim:=proc(r1,r2,r12,n,D1,E,Z,MaxD) local gu,F,R,X,x,y,i:

X:=[seq(x[i],i=1..n),seq(y[i],i=1..n)]:
F:=[sqrt(add(x[i]^2,i=1..n)),sqrt(add(y[i]^2,i=1..n)),
sqrt(add((x[i]-y[i])^2,i=1..n))]:
R:=[r1,r2,r12]:

gu:=Chaim1(x,y,n,D1,E,Z):
gu:=ConOp(F,R,X,D1,gu,MaxD):


gu:
end:



#GuessPol(f,X,F,R,MaxD): Given an expression f in the variables
#X[1], ..., X[n] that is known to be a polynomial in
#R[1],...,R[m] with R[i] substituted for F[i] where
#the F[i] is an expression in the X[1],...,X[n]
#tries to find that polynomial in R[1],..., R[m]
#For example, try:
#GuessPol(sqrt(x^2+y^2)^3,[x,y],[sqrt(x^2+y^2)],[r],5);
GuessPol:=proc(f,X,F,R,MaxD) local d1,gu:

for d1 from 0 to MaxD do
 gu:=GuessPol1(f,X,F,R,d1):
 
 if gu<>FAIL then
  RETURN(gu):
 fi:

od:

FAIL:
end:


#Comps(d,k): all the vectors of length k that add-up to d exactly
#For example, try:
#Comps(3,2);
Comps:=proc(d,k) local gu1,gu,i,g1:
option remember:

if k=0 then
 if d=0 then
  RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:

for i from 0 to d do
 gu1:=Comps(d-i,k-1):
 gu:=gu union {seq([op(g1),i],g1 in gu1)}:
od:
gu:
end:

  


#CompsA(d,k): all the vectors of length k that add-up to <=d 
#For example, try:
#CompsA(3,2);
CompsA:=proc(d,k) local d1: {seq(op(Comps(d1,k)),d1=0..d)}: end:







#GenPol(X,d,a): a generic polynomial in the list of variables X
#of degree d with coefficients indexed by the symbol a,
#followed by the set of coefficients
#For example, try:
#GenPol([x,y],3,a);
GenPol:=proc(X,d,a) local gu,g,k,i:
k:=nops(X):
gu:=CompsA(d,k):
add(a[op(g)]*mul(X[i]^g[i],i=1..k),g in gu),
{seq(a[op(g)],g in gu)}:
end:


#GuessPol1(f,X,F,R,d1): Given an expression f in the variables
#X[1], ..., X[n] that is known to be a polynomial in
#R[1],...,R[m] with R[i] substituted for F[i] where
#the F[i] is an expression in the X[1],...,X[n]
#tries to find that polynomial in R[1],..., R[m] of degree<=d1 .
#For example, try:
#GuessPol1(sqrt(x^2+y^2)^3,[x,y],[sqrt(x^2+y^2)],[r],3);
GuessPol1:=proc(f,X,F,R,d1) local pol,a,var,lu,f1,i1,L,eq:
pol:=GenPol(R,d1,a):
var:=pol[2]:
pol:=pol[1]:
eq:={}:


lu:=CompsA(d1+2,nops(X)):
if nops(lu)-nops(var)<5 then
 print(`not enough equations`):
 RETURN(FAIL):
fi:

f1:=f-subs({seq(R[i1]=F[i1],i1=1..nops(R))},pol):

eq:={seq(subs({seq(X[i1]=L[i1],i1=1..nops(X))},f1),L in lu)}:


var:=solve(eq,var):

if var=NULL then
FAIL:
else
subs(var,pol):
fi:


end:


#GuessPol(f,X,F,R,MaxD): Given an expression f in the variables
#X[1], ..., X[n] that is known to be a polynomial in
#R[1],...,R[m] with R[i] substituted for F[i] where
#the F[i] is an expression in the X[1],...,X[n]
#tries to find that polynomial in R[1],..., R[m]
#For example, try:
#GuessPol(sqrt(x^2+y^2)^3,[x,y],[sqrt(x^2+y^2)],[r],5);
GuessPol:=proc(f,X,F,R,MaxD) local d1,gu:

for d1 from 0 to MaxD do
 gu:=GuessPol1(f,X,F,R,d1):
 
 if gu<>FAIL then
  RETURN(gu):
 fi:

od:

FAIL:
end:

#YafePlus(P,Y,X,F,R,MaxD): given a multivariable polynomial P
#in Y[1],..,Y[k] (Y is a list of variables)
#collects and simplifies
#where the coeffs. are functions of X[1], ..., X[n]
#believed to be expressible as rational functions
#in the R[1], .., R[k] where R[1]=F[1], ...
#for example, try
#YafePlus(x^2*z+2*x*y*z+y^2*z,[z],[x,y],[x+y],[r]);
YafePlus:=proc(P,Y,X,F,R,MaxD) local i,lu,P1,v,gu,i1,gu1,L,gu11,gu12,
mu11,mu12:

if not type(P,`+`) then
 RETURN(normal(P)):
fi:

lu:={}:

for i from 1 to nops(P) do
P1:=op(i,P):
v:=[seq(degree(P1,Y[i1]),i1=1..nops(Y))]:
lu:=lu union {v}:
od:

gu:=0:

for L in lu do
gu1:=P:

for i1 from 1 to nops(Y) do
 gu1:=coeff(gu1,Y[i1],L[i1]):
od:

gu11:=numer(gu1):
gu12:=denom(gu1):


mu11:=GuessPol(gu11,X,F,R,MaxD):


if mu11<>FAIL then
 gu11:=mu11:
fi:

mu12:=GuessPol(gu12,X,F,R,MaxD):

if mu12<>FAIL then
 gu12:=mu12:
fi:


gu:=gu+ factor(normal(gu11/gu12))*mul(Y[i1]^L[i1],i1=1..nops(Y)):
od:

gu:

end:

#ConOp(F,R,X,D1,Ope,MaxD): 
#Given a transormation 
#X[1]=F[1](R[1],..R[m]), X[2]=F[2](R[1],..R[m])
#expressing n X's in terms of m R's, where R is a list
#of independent variables and X is a list of dependent variables
#and D1 is a letter indicating that we denote diff. w.r.t. a variable
#v by D1[v], and Ope is a linear differential operator
#in the X's and D1[X]'s converts it to a diff. operator
#in D1[R[1]], ..., D1[R1[m]] with coeff. that are 
#functions of X[1], ..., X[n]
#For example, try:
#ConOp([x^2+y^2,x*y],[u,v],[x,y],E,E[x]^2+E[y]^2,4);
ConOp:=proc(F,R,X,D1,Ope,MaxD) local i,gu:

if not type(Ope,`+`) then
 RETURN(ApplyMono(F,R,X,D1,Ope)):
fi:

gu:=expand(add(ApplyMono(F,R,X,D1,op(i,Ope)),i=1..nops(Ope))):


YafePlus(gu,[seq(D1[R[i]],i=1..nops(R))],X,F,R,MaxD):

end:





#Hafukh(F,X,R): Given an affine linear transformation
#from a list of variables X to a list of variables R
#given by F, returns the inverse transformation
#For example, try:
#Hafukh([x+3*y,2*x-y],[x,y],[u,v]);

Hafukh:=proc(F,X,R) local i,var:
if nops(F)<>nops(R) then
 print(`F and R should have the same length`):
 RETURN(FAIL):
fi:

var:=solve({seq(F[i]=R[i],i=1..nops(F))},convert(X,set)):

if var=NULL then
 RETURN(FAIL):
fi:

[seq(subs(var,X[i]),i=1..nops(X))]:
end:



#ConOpLin(F,R,X,D1,Ope): 
#Given a linear transormation 
#X[1]=F[1](R[1],..R[m]), X[2]=F[2](R[1],..R[m])
#expressing n X's in terms of m R's, where R is a list
#of independent variables and X is a list of dependent variables
#and D1 is a letter indicating that we denote diff. w.r.t. a variable
#v by D1[v], and Ope is a linear differential operator
#in the X's and D1[X]'s converts it to a diff. operator
#in D1[R[1]], ..., D1[R1[m]] with coeff. that are 
#functions of X[1], ..., X[n],
#convers that operator into a linear differential operator
#in the variables R[1], ..., R[m]
#For example, try:
#ConOpLin([x+3*y,2*x-y],[u,v],[x,y],E,E[x]^2+E[y]^2,4);
ConOpLin:=proc(F,R,X,D1,Ope) local i,gu,lu:



if not type(Ope,`+`) then
gu:=ApplyMono(F,R,X,D1,Ope):
lu:=Hafukh(F,X,R):
gu:=subs({seq(X[i]=lu[i],i=1..nops(X))},gu):
RETURN(gu):
fi:



gu:=expand(add(ApplyMono(F,R,X,D1,op(i,Ope)),i=1..nops(Ope))):

lu:=Hafukh(F,X,R):
gu:=subs({seq(X[i]=lu[i],i=1..nops(X))},gu):
gu:=Yafe(gu,[seq(D1[R[i]],i=1..nops(R))]):


end:


ChaimPC3:=proc(r1,r2,r12,D1,E,Z):
2*(Z*r2*r12+Z*r1*r12-r1*r2+E*r1*r2*r12)/r1/r2/r12+4/r12*D1[r12]+2*D1[r12]^2+2/
r2*D1[r2]-(-r12^2-r2^2+r1^2)/r2/r12*D1[r2]*D1[r12]+D1[r2]^2+2/r1*D1[r1]+(r12^2-
r2^2+r1^2)/r1/r12*D1[r1]*D1[r12]+D1[r1]^2:
end:

ChaimPC:=proc(r1,r2,r12,n,D1,E,Z):
2*(Z*r2*r12+Z*r1*r12-r1*r2+E*r1*r2*r12)/r1/r2/r12+4/r12*D1[r12]+2*D1[r12]^2+2/
r2*D1[r2]-(-r12^2-r2^2+r1^2)/r2/r12*D1[r2]*D1[r12]+D1[r2]^2+2/r1*D1[r1]+(r12^2-
r2^2+r1^2)/r1/r12*D1[r1]*D1[r12]+D1[r1]^2
+(n-3)*(2*D1[r12]/r12+D1[r1]/r1+D1[r2]/r2):
:
end:


#Pekeris1(u,v,w,D1,e,Z,n): The diff. operator Chaim in terms of the perimetric
#coordinates u,v,w, and D1 denotes differentiation in n dimensions. Do:
#Pekeris1(u,v,w,D1,e,Z,n);
Pekeris1:=proc(u,v,w,D1,e,Z,n) local gu,E,r1,r2,r12:
gu:=ChaimPC(r1,r2,r12,n,D1,E,Z):
gu:=expand(subs(E=-e^2,gu)/4/e):
gu:=ConOpLin([e*(r2+r12-r1),e*(r1+r12-r2),2*e*(r1+r2-r12)],
[u,v,w],[r1,r2,r12],D1,gu):
gu:=expand(gu*(2*v+w)*(w+2*u)*(u+v)):
gu:=subs({D1[u]=D1[u]-1/2,D1[v]=D1[v]-1/2,D1[w]=D1[w]-1/2},gu):
gu:=normal(expand(gu)):
RETURN(gu):
gu:=Yafe(gu,[D1[u],D1[v],D1[w]]):

gu:
end:


#Kefel(ope1,ope2,n,N): multiplies ope1(n,N) by ope2(n,N)
#For example, try:
#Kefel(N,n,n,N);
Kefel:=proc(ope1,ope2,n,N) local i,gu:
gu:=0:

for i from ldegree(ope1,N) to degree(ope1,N) do
 gu:=gu+coeff(ope1,N,i)*subs(n=n+i,ope2)*N^i:
od:
expand(gu):
end:


#Lag(n,x):The Laguerre polynomials L_n(x)
Lag:=proc(n,x) local k:
add(binomial(n,k)*(-x)^k/k!,k=0..n):
end:


#SimplifyMonomial1(L,x,D,n,N): Given a list L of
#Rewriting Rules, where L[i] is an expression
#for x^(i-1)*D^*(i-1) as an operator in n, N, 
#(up to nops(L))
#finds the simplified form an operator in n,N
#only for the Monomial of the x^a*D^b with a>=b>=nops(L)-1
#For example, try:
#SimplifyMonomial1([-(n+1)*N+(2*n+1)-n/N,n-n/N,(n^2-n)*(N-2/N+1/N^2)],x,D,n,N,x^3*D^2);
SimplifyMonomial1:=proc(L,x,D,n,N,Mono) local a,b,ope,k,ope1,fu:
a:=degree(Mono,x):
b:=degree(Mono,D):

if a<b or b>=nops(L) then
 print(`bad input`):
 RETURN(FAIL):
fi:

fu:=normal(Mono/x^a/D^b):

if diff(fu,x)<>0 or diff(fu,D)<>0 then
 RETURN(FAIL):
fi:

if a=0 and b=0 then
 RETURN(fu):
fi:

if b=0 then
 ope1:=L[1]:
 ope:=L[1]:
 
 for k from 1 to a-1 do
   ope:=Kefel(ope1,ope,n,N):
  od:
RETURN(fu*ope):
fi:


ope:=L[b+1]:
ope1:=L[1]:

for k from 1 to a-b do
ope:=Kefel(ope1,ope,n,N):
od:
expand(ope):
end:







KlmnWithTypo:=proc(l,m,n) local om:
om:=l+m+n:
(1/24)*om*(om+2)*(2*om+5)+1/16*(1-(-1)^om)+
(l+m)^2/4+(1-(-1)^(l+m))/8+l+1:
end:

#The linearization of the vectors (l,m,n) with l<=m
#For example, try: Klmn(0,1,1);
Klmn:=proc(l,m,n) local om:
om:=l+m+n:
(1/24)*om*(om+2)*(2*om+5)+1/16*(1-(-1)^om)+
(l+m)^2/4+(l+m)/2+(1-(-1)^(l+m))/8+l+1:
end:





#ConvertMonomial(u,v,w,l,m,n,L,M,N,Li,Mono): Given
#a monomial in u,v,w, Mono, and discrete variables
#l,m,n, corresponding to them(resp.) and a symbol S1
#for the shift and a list Li, of length 3, of operators representating
#left multiplication by u,v,w resp. finds the corresponding
#recurrence operator in terms of L,M,N
#For example try:
#ConvertMonomial(u,v,w,l,m,n,S1,[-(l+1)*L+2*l+1-l*L^(-1),-(m+1)*M+2*m+1-m*M^(-1),-(n+1)*N+2*n+1-n*N^(-1)],u*v*w):
ConvertMonomial:=proc(u,v,w,l,m,n,L,M,N,Li,Mono) local a,b,c,lu,
gu,i,A,B,C,lu1,c1,c2,c3:

a:=degree(Mono,u):
b:=degree(Mono,v):
c:=degree(Mono,w):

lu:=normal(Mono/u^a/v^b/w^c):

A:=degree(lu,l):
B:=degree(lu,m):
C:=degree(lu,n):


lu1:=normal(lu/l^A/m^B/n^C):


if {normal(diff(lu1,u)),normal(diff(lu1,w)),normal(diff(lu1,v))}<>{0} then
  RETURN(FAIL):
fi:

c1:=degree(lu1,L):
c2:=degree(lu1,M):
c3:=degree(lu1,N):

lu1:=lu1/L^c1/M^c2/N^c3:
gu:=lu1:

for i from 1 to a do
gu:=Kefel(Li[1],gu,l,L):
od:

for i from 1 to b do
gu:=Kefel(Li[2],gu,m,M):
od:

for i from 1 to c do
gu:=Kefel(Li[3],gu,n,N):
od:


gu:=l^A*m^B*n^C*subs({l=l+c1,m=m+c2,n=n+c3},gu)*L^c1*M^c2*N^c3:
expand(gu):

end:

#ConvertOp(u,v,w,l,m,n,L,M,N,Li,Ope): Given
#an operator in u,v,w, Ope, and discrete variables
#l,m,n, corresponding to them(resp.) and a symbol S1
#for the shift and a list Li, of length 3, of operators representating
#left multiplication by u,v,w resp. finds the corresponding
#recurrence operator in terms of S1[l],S1[m],S1[n]
#For example try:
#ConvertOp(u,v,w,l,m,n,S1,[-(l+1)*L+2*l+1-l*L^(-1),-(m+1)*M+2*m+1-m*M^(-1),-(n+1)*N+2*n+1-n*N^(-1)],u+v+w);
ConvertOp:=proc(u,v,w,l,m,n,L,M,N,Li,Ope) local i:

if not type(Ope,`+`) then
RETURN(ConvertMonomial(u,v,w,l,m,n,L,M,N,Li,Ope)):
fi:

expand(add(ConvertMonomial(u,v,w,l,m,n,L,M,N,Li,op(i,Ope)),i=1..nops(Ope))):
end:


#Pekeris2(l,m,n,e,Z,L,M,N,n): The recurrence
#operator corresponding to Pekeris1
#and using shift operators corresponding to u,v,w
#in n dimensions
#Do:
#Pekeris2(l,m,n,e,Z,L,M,N,3);
Pekeris2:=proc(l,m,n,e,Z,L,M,N,memad) local gu,u,v,w,D1:
option remember:
gu:=Pekeris1(u,v,w,D1,e,Z,memad):
gu:=subs(
{
D1[u]^2=(u-1)/u*D1[u]-l/u,
D1[v]^2=(v-1)/v*D1[v]-m/v,
D1[w]^2=(w-1)/w*D1[w]-n/w
},gu);
gu:=expand(gu):
gu:=subs({
D1[u]=l/u-l/u/L,
D1[v]=m/v-m/v/M,
D1[w]=n/w-n/w/N
},gu):
gu:=expand(gu):
gu:=
ConvertOp(u,v,w,l,m,n,L,M,N,
[-(l+1)*L+2*l+1-l*L^(-1),
-(m+1)*M+2*m+1-m*M^(-1),
-(n+1)*N+2*n+1-n*N^(-1)],gu);
gu:
gu:=Yafe(gu,[L,M,N]):
end:



#Pek1h(l,m,n,e,Z,L,M,N): The recurrence
#operator corresponding to Pek1
#and using shift operators corresponding to u,v,w
#Do:
#Pek1h(l,m,n,e,Z,L,M,N);
Pek1h:=proc(l,m,n,e,Z,L,M,N) local gu,u,v,w,D1:
gu:=Pek1(u,v,w,D1,e,Z):
gu:=subs(
{
D1[u]^2=(u-1)/u*D1[u]-l/u,
D1[v]^2=(v-1)/v*D1[v]-m/v,
D1[w]^2=(w-1)/w*D1[w]-n/w
},gu);
gu:=expand(gu):
gu:=subs({
D1[u]=l/u-l/u/L,
D1[v]=m/v-m/v/M,
D1[w]=n/w-n/w/N
},gu):
gu:=expand(gu):
RETURN(gu):
gu:=Yafe(gu,[L,M,N]):
end:


#MakeEqs(A,Z,e,om,memad): makes the set of equations for
#Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om
#in memad dimensions
#The output is the set of equations followed by the
#set of variables.For example, do:
#MakeEqs(A,Z,e,1,3);
MakeEqs:=proc(A,Z,e,om,memad) local eq,var,ope,L,M,N,l,m,n,l1,m1,n1,
i1,i2,i3,eq1:
option remember:
ope:=Pekeris2(l,m,n,e,Z,L,M,N,memad):
var:={}:
eq:={}:

for l1 from 0 to om do
 for m1 from 0 to om-l1 do
  for n1 from 0 to om-l1-m1 do
  eq1:=0:
   for i1 from ldegree(ope,L) to degree(ope,L) do
   for i2 from ldegree(ope,M) to degree(ope,M) do
   for i3 from ldegree(ope,N) to degree(ope,N) do
   if l1+i1>=0 and  m1+i2>=0 and  n1+i3>=0 and
    (l1+i1)+(m1+i2)+(n1+i3)<=om then
  eq1:=eq1+
  subs({l=l1,m=m1,n=n1},coeff(coeff(coeff(ope,L,i1),M,i2),N,i3))*
     A[l1+i1,m1+i2,n1+i3]:
   fi:
   od:
  od:
od:
eq:=eq union {eq1}:
var:=var union {A[l1,m1,n1]}:
od:
od:
od:

eq,var:
end:



#MakeEqsPara(A,Z,e,om,memad): makes the set of equations for
#Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om,
#with the para assumption (i.e. A(l,m,n)=A(m,l,n))
#The output is the set of equations followed by the
#set of variables.For example, do:
#MakeEqsPara(A,Z,e,1,3);
MakeEqsPara:=proc(A,Z,e,om,memad) local gu,l1,m1,n1,eq,var:
option remember:
gu:=MakeEqs(A,Z,e,om,memad):
eq:=gu[1]:
var:=gu[2]:

for n1 from 0 to om do
for l1 from 0 to om-n1 do
for m1 from 0 to om-n1-l1 do
 if l1>m1 then
   eq:=subs(A[l1,m1,n1]=A[m1,l1,n1],eq):
   var:=var minus {A[l1,m1,n1]}:
 fi:
od:
od:
od:
eq,var:

end:



#MatrixPara(Z,e,om,memad): makes the matrix of
#Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om,
#with the para assumption (i.e. A(l,m,n)=A(m,l,n))
#in memad dimensions
#The output is the set of equations followed by the
#set of variables.For example, do:
#MatrixPara(Z,e,5,3);
MatrixPara:=proc(Z,e,om,memad) local gu,eq,var,i,j,A:
option remember:
gu:=MakeEqsPara(A,Z,e,om,memad) :
eq:=gu[1]:
var:=gu[2]:
var:=convert(var,list):
[seq([seq(coeff(eq[i],var[j]),j=1..nops(var))],i=1..nops(eq))]:
end:

with(linalg):

#FastDet(lu,e): computes a symbolic determint of
# a matrix lu with entries that depend on variable e
#For example, try:
#FastDet([[e,1],[1,e]],e);
FastDetSlow:=proc(lu,e) local gu,pol,a,i,eq,var:
gu:=[seq(det(subs(e=i,lu)),i=1..nops(lu)+2)]:
pol:=add(a[i]*e^i,i=0..nops(lu)):
var:={seq(a[i],i=0..nops(lu))}:
eq:={seq(subs(e=i,pol)-gu[i],i=1..nops(lu)+2)}:
var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:
sort(subs(var,pol)):
end:







FastDet1 := proc(A)
local var, deg, det, det1, n, i, j, p, pp, vals, a, b, len;
   n := {LinearAlgebra[Dimension](A)};
   if nops(n) <> 1 then return FAIL fi:
#"Matrix is not square." fi;
   n := op(1,n);
   var := indets(A);
   if nops(var) <> 1 then return FAIL: fi:
#"Matrix entries must be univariate polynomials." fi;
   var := op(1,var);
   deg := add(max(seq(degree(A[i,j]),i=1..n)),j=1..n);
   a := ceil(-deg/2); b := ceil(deg/2);
   det := 0; det1 := 1;
   p := 2^24; pp := 1;
   while det <> det1 do
     det1 := det;
     p := prevprime(p);
     vals := [seq(Det(subs(var=i,A)) mod p, i=a..b)];
     det := Interp([seq(i,i=a..b)], vals, var) mod p;
     det := [seq(coeff(det, var, i), i=0..deg)];
     if det1 <> 1 then
       len := max(nops(det), nops(det1));
       det := [op(det), seq(0, i=1..len-nops(det))];
       det1 := [op(det1), seq(0, i=1..len-nops(det1))];
       det := [seq(chrem([det[i],det1[i]], [p,pp]), i=1..len)];
     fi;
     pp := p*pp;
     det := mods(det, pp);
   od;
   return add(det[i+1]*var^i, i=0..deg);
end:





#MatrixOrtho(Z,e,om,memad): makes the matrix of
#Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om,
#with the ortho assumption (i.e. A(l,m,n)=-A(m,l,n))
#The output is the set of equations followed by the
#set of variables.For example, do:
#MatrixOrtho(Z,e,1,3);
MatrixOrtho:=proc(Z,e,om,memad) local gu,eq,var,i,j,A:
option remember:
gu:=MakeEqsOrtho(A,Z,e,om,memad) :
eq:=gu[1]:
var:=gu[2]:
var:=convert(var,list):
[seq([seq(coeff(eq[i],var[j]),j=1..nops(var))],i=1..nops(eq))]:
end:


#MakeEqsOrtho(A,Z,e,om,memad): In memad dimensions
#makes the set of equations for
#Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om,
#with the ortho assumption (i.e. A(l,m,n)=-A(m,l,n))
#The output is the set of equations followed by the
#set of variables.For example, do:
#MakeEqsOrtho(A,Z,e,1);
MakeEqsOrtho:=proc(A,Z,e,om,memad) local gu,l1,m1,n1,eq,var,eqN:
option remember:
gu:=MakeEqs(A,Z,e,om,memad):
eq:=gu[1]:
var:=gu[2]:

for n1 from 0 to om do
for l1 from 0 to om-n1 do
for m1 from 0 to om-n1-l1 do
 if l1>m1 then
   eq:=subs(A[l1,m1,n1]=-A[m1,l1,n1],eq):
   var:=var minus {A[l1,m1,n1]}:
elif l1=m1 then
   eq:=subs(A[l1,m1,n1]=0,eq):
   var:=var minus {A[l1,m1,n1]}:
 fi:
od:
od:
od:


eq:=eq minus {0}:
eqN:={}:

while eq<>{} do
eqN:=eqN union {eq[1]}:
eq:=eq minus {eq[1],-eq[1]}:
od:
eqN,var:
end:



FastDet := proc(A)
FastDet1(Matrix(A)):
end:

#FindRoot1(A,e,L): inputs a matrix A whose entries are expressions
#and e and that has an eigenvalue in the interval L (of the form
#[a,b]), outputs the interval [a,(a+b)/2] or [(a+b),b] that
#has it. Try:
#FindRoot1(MatrixPara(2,e,5,3),e,[17/10,18/10]);
FindRoot1:=proc(A,e,L) local lu1,lu2,lu,c,a,b:
a:=L[1]:
b:=L[2]:
lu1:=det(subs(e=a,A)):
lu2:=det(subs(e=b,A)):

 if lu1*lu2>0 then
   RETURN(FAIL):
 fi:
c:=(a+b)/2:
lu:=det(subs(e=c,A)):

if lu1*lu<0 then
 RETURN([a,c]):
elif lu=0 then
  RETURN([c,c]):
else
  RETURN([c,b]):
fi:

end:




#FindRoot(A,e,L,prec): inputs a matrix A whose entries are expressions
#and e and that has an eigenvalue in the interval L (of the form
#[a,b]), and the required precision
#outputs the interval [a,(a+b)/2] or [(a+b),b] that
#has it. Try:
#FindRoot(MatrixPara(2,e,5,3),e,[17/10,18/10],10);
FindRoot:=proc(A,e,L,prec) local L1:

L1:=L:

while L1[2]-L1[1]>10^(-prec) do
 L1:=FindRoot1(A,e,L1,prec):
od:

L1:
end:




EVold:=proc(C,e) local A,B,i,j:

A:=[seq([seq(coeff(C[i][j],e,0),j=1..nops(C[i])) ],i=1..nops(C))]:
B:=[seq([seq(coeff(C[i][j],e,1),j=1..nops(C[i])) ],i=1..nops(C))]:
A:=Matrix(A):
B:=Matrix(B):
LinearAlgebra[Eigenvalues](A.LinearAlgebra[MatrixInverse](B)):
end:


#EV(C,e): Given a square-matrix C whose entries are affine-linear
#expressions in e, finds all the eigenvalues
#For example, try:
#EV(MatrixPara(2,e,5,3));
EV:=proc(C,e) local A,B,i,j,T:
option remember:
A:=[seq([seq(evalf(coeff(C[i][j],e,0)),j=1..nops(C[i])) ],i=1..nops(C))]:
B:=[seq([seq(evalf(coeff(C[i][j],e,1)),j=1..nops(C[i])) ],i=1..nops(C))]:
A:=Matrix(A):
B:=Matrix(B):
T:=LinearAlgebra[Eigenvalues](A,B):
[seq(-coeff(T[i],I,0),i=1..nops(C))]:
end:



sn:=proc(resh,n):
-1/log(op(n+1,resh)*op(n-1,resh)/op(n,resh)^2):
end:

#Zinn(resh,n)'s method: For example, try:
#Zinn([seq(1.1^i*(i+1),i=1..20)],19);
Zinn:=proc(resh,n)
local s1,s2:
s1:=sn(resh,n):
s2:=sn(resh,n-1):
evalf(2*(s1+s2)/(s1-s2)^2),
evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
end:
 
 
 

#InterpolateCP(L): given a list L of 4 consecutive values 
#guesses the limiting value using eqs. (33) and (34) in
#Pekeris' paper
InterpolateCP:=proc(L) local x,a,mu:

mu:=
solve({a*x=(L[3]-L[2])/(L[2]-L[1]), (a+1)/2*x=(L[4]-L[3])/(L[3]-L[2])},
{a,x}):

if mu=NULL then
 RETURN(L[4]):
fi:

a:=subs(mu,a):
x:=subs(mu,x):

coeff(L[1]+(L[2]-L[1])/(1-x)^a,I,0):

end:



#EnergyPara(omega,Z,d,memad,orekh): 
EnergyPara:=proc(omega,Z0,d,memad,orekh) local gu,i,L,cha,e,C,Z:
if omega-(orekh-1)*d<=2 or orekh<=4 then
 print(`omega-(orekh-1)*d must be >=3 and orekh>=5 `):
 RETURN(FAIL):
fi:

L:=[seq(EVf(subs(Z=Z0,MatrixPara(Z,e,omega+i*d,memad)),e),i=-(orekh-1)..0)]:
cha:=[seq(InterpolateCP([op(j..j+3,L)]),j=1..nops(L)-3)]:
cha[nops(cha)],L, cha, cha[nops(cha)]-cha[nops(cha)-1],L[nops(L)]-L[nops(L)-1]:
end:


#EnergyOrtho(omega,Z,d,memad,orekh): 
EnergyOrtho:=proc(omega,Z,d,memad,orekh) local gu,i,L,cha,e:
if omega-(orekh-1)*d<=2 or orekh<=4 then
 print(`omega-(orekh-1)*d must be >=3 and orekh>=5 `):
 RETURN(FAIL):
fi:

L:=[seq(EVf(MatrixOrtho(Z,e,omega+i*d,memad),e),i=-(orekh-1)..0)]:
cha:=[seq(InterpolateCP([op(j..j+3,L)]),j=1..nops(L)-3)]:
cha[nops(cha)],L, cha, cha[nops(cha)]-cha[nops(cha)-1],L[nops(L)]-L[nops(L)-1]:
end:





#Adj11(ope,n,N): the adjoint of the operator ope(n,N)
Adj11:=proc(ope,n,N) local i:
add(factor(subs(n=n-i,coeff(ope,N,i)))*N^(-i),i=ldegree(ope,N)..degree(ope,N)):
end:

Adj1:=proc(ope,Lin,LiN) local i,gu:

if nops(Lin)<>nops(LiN) then
RETURN(FAIL):
fi:

gu:=ope:

for i from 1 to nops(Lin) do
 gu:=Adj11(gu,Lin[i],LiN[i]):
od:

gu:=expand(gu):

Yafe(gu,LiN):

end:


#EVf(C,e): Given a square-matrix C whose entries are affine-linear
#expressions in e, finds all the eigenvalues
#For example, try:
#EVf(MatrixPara(2,e,5,3));
EVf:=proc(C,e) local A,B,i,j,T,C1:
option remember:
A:=[seq([seq(evalf(-coeff(C[i][j],e,0)),j=1..nops(C[i])) ],i=1..nops(C))]:
B:=[seq([seq(evalf(-coeff(C[i][j],e,1)),j=1..nops(C[i])) ],i=1..nops(C))]:
A:=Matrix(A):
B:=Matrix(B):
C1:=LinearAlgebra[Multiply](LinearAlgebra[MatrixInverse](B),A):

T:=LinearAlgebra[Eigenvalues](C1):
-coeff(T[1],I,0):
end:






#PaperPara(psi,x,y,r,u,v,w,F,MaxZ,omega,A,E,e,Z,m,l,n,d,orekh):
#outputs a paper about the ground (and excited) energy
#psi,x,y,r,u,v,w,F,A,E,e,Z,m,l,n are symbols as in Pekeris'
#paper, MaxZ is the the maximal charge (i.e. Z is taken from 1 to MaxZ)
# omega is the largest omega considered
#For example, to reproduce
#Pekeris' paper type:
#     PaperPara(psi,x,y,r,u,v,w,F,10,11,A,E,e,Z,l,m,n,1,4);`):
PaperPara:=proc(psi,x,y,r,u,v,w,F,MaxZ,omega,A,E,e,Z,l,m,n,d,orekh)
local D,ope1,L,ope2,lu,i1,i2,i3,Z0,vu,t0,t1:
t0:=time():

print(`Ground State of Two-Electron Atoms (Para states) according to Pekeris`):
print(``):
print(`Recall that the two-electron Scroedinger Equation for`):
print(`the state function`, psi(x[1],x[2],x[3],y[1],y[2],y[3])):
print(`in atomic units, satisfies the Schroedinger equation`):
print((Chaim1a(x,y,3,D,E,Z))*psi=0):
print(`where r1,r2 are the distances of electron 1 and 2 from the origin`):
print(`and r12 is the distance between them `):
print(``):
print(`Since the ground state only depends on r1,r2,r12`):
print(`We see that, by changing coordinates that`, psi(r1,r2,r12)):
print(`satisfies the differential equation`):
print(ChaimPC(r1,r2,r12,3,D,E,Z)*psi(r1,r2,r12)=0):
print(`Let E=-e^2 and `):
print(`Introducing so-called perimetric coordinates`):
print(u=e*(r2+r12-r1),v=e*(r1+r12-r2),w=2*e*(r1+r2-r12)):
print(`and making the change of dependent variable`):
print(psi=exp(-(u+v+w)/2)*F(u,v,w)):
print(`we see that`, F(u,v,w),  `is annihilated by`):
print(`the linear partial diferential operator:`):
ope1:=Pekeris1(u,v,w,D,e,Z,3):
print(collect(ope1,[D[u],D[v],D[w]])):
print(`writing`, F(u,v,w), `in terms of Laguerre polynomials`):
print(F(u,v,w)=Sum(Sum(Sum(A[l,m,n]*L[l](u)*L[m](v)*L[k](w),l=0..infinity),m=0..infinty),n=0..infinity)):
print(`and using the well-known relations (easily provable by WZ)`):
print(x*diff(L[n](x),x$2)=(x-1)*diff(L[n](x),x)-n*L[n](x)):
print(``):
print(x*diff(L[n](x),x)=n*L[n](x)-(n-1)*L[n-1](x)):
print(``):
print(x*L[n](x)=-(n+1)*L[n+1](x)+(2*n+1)*L[n](x)-n*L[n-1](x)):
print(``):
ope2:=Pekeris2(l,m,n,e,Z,L,M,N,3):
print(`We get that the coefficients`, A[l,m,n] , `are annihilated`):
print(`by the linear partial recurrence equation with `):
print(`polynomial coefficients`):

print(collect(ope2,[L,M,N])):
print(` where `, L,M,N, `are the shift operators in the discrete variables`):
print(l,m,n, `respectively. `):

lu:=0:

for i1 from ldegree(ope2,L) to degree(ope2,L) do
for i2 from ldegree(ope2,M) to degree(ope2,M) do
for i3 from ldegree(ope2,N) to degree(ope2,N) do
lu:=lu+
coeff(coeff(coeff(ope2,L,i1),M,i2),N,i3)*A[l+i1,m+i2,n+i3]:
od:
od:
od:
print(`Or more explicitly, in humanese`):
print(lu=0):

t1:=time():
print(`This ends the symbolic part, meticulously done by C.L.Pekeris`):
print(`It took`, t1-t0, `seconds `):
print(``):
print(`Now is time for the numeric part, done by WEIZAC and`):
print(`cleverly programmed in machine language by Yigal Accad,`):
print(`taking account that we are talking about para states`):
print(`so that A[l,m,n]=A[m,l,n], `):

print(``):
print(`By truncating this infinite system to those`, l,m,n):
print(`satisfying `, l+m+n<=omega , `we get that the parameter`, e):
print(`where E=-e^2 ,for charges Z between 1, and `, MaxZ, `are as follows`):



for Z0 from 1 to MaxZ do
vu:=EnergyPara(omega,Z0,d,3,orekh):
print(`When the charge`,Z=Z0, `the appx. value is`,e=vu[1]):
print(`For omega=`,[seq(omega+i*d,i=-(orekh-1)..0)], `the e's for the truncated matrix are, resp.`):
print(vu[2]):
print(`The extrapolated values are`, vu[3]):
print(`the probable error is`, vu[4]):
od:

print(``):
print(`The whole thing took`, time()-t0, `seconds `):
end:



#PaperOrtho(psi,x,y,r,u,v,w,F,MaxZ,omega,A,E,e,Z,m,l,n,d,orekh):
#outputs a paper about the ground (and excited) energy
#psi,x,y,r,u,v,w,F,A,E,e,Z,m,l,n are symbols as in Pekeris'
#paper, MaxZ is the the maximal charge (i.e. Z is taken from 1 to MaxZ)
# omega is the largest omega considered
#For example, to reproduce
#Pekeris' paper type:
#     PaperOrtho(psi,x,y,r,u,v,w,F,10,11,A,E,e,Z,l,m,n,1,4);`):
PaperOrtho:=proc(psi,x,y,r,u,v,w,F,MaxZ,omega,A,E,e,Z,l,m,n,d,orekh)
local D,ope1,L,ope2,lu,i1,i2,i3,Z0,vu,t0,t1:
t0:=time():

print(`Ground State of Two-Electron Atoms (Ortho states) according to Pekeris`):
print(``):
print(`Recall that the two-electron Scroedinger Equation for`):
print(`the state function`, psi(x[1],x[2],x[3],y[1],y[2],y[3])):
print(`in atomic units, satisfies the Schroedinger equation`):
print((Chaim1a(x,y,3,D,E,Z))*psi=0):
print(`where r1,r2 are the distances of electron 1 and 2 from the origin`):
print(`and r12 is the distance between them `):
print(``):
print(`Since the ground state only depends on r1,r2,r12`):
print(`We see that, by changing coordinates that`, psi(r1,r2,r12)):
print(`satisfies the differential equation`):
print(ChaimPC(r1,r2,r12,3,D,E,Z)*psi(r1,r2,r12)=0):
print(`Let E=-e^2 and `):
print(`Introducing so-called perimetric coordinates`):
print(u=e*(r2+r12-r1),v=e*(r1+r12-r2),w=2*e*(r1+r2-r12)):
print(`and making the change of dependent variable`):
print(psi=exp(-(u+v+w)/2)*F(u,v,w)):
print(`we see that`, F(u,v,w),  `is annihilated by`):
print(`the linear partial diferential operator:`):
ope1:=Pekeris1(u,v,w,D,e,Z,3):
print(collect(ope1,[D[u],D[v],D[w]])):
print(`writing`, F(u,v,w), `in terms of Laguerre polynomials`):
print(F(u,v,w)=Sum(Sum(Sum(A[l,m,n]*L[l](u)*L[m](v)*L[k](w),l=0..infinity),m=0..infinty),n=0..infinity)):
print(`and using the well-known relations (easily provable by WZ)`):
print(x*diff(L[n](x),x$2)=(x-1)*diff(L[n](x),x)-n*L[n](x)):
print(``):
print(x*diff(L[n](x),x)=n*L[n](x)-(n-1)*L[n-1](x)):
print(``):
print(x*L[n](x)=-(n+1)*L[n+1](x)+(2*n+1)*L[n](x)-n*L[n-1](x)):
print(``):
ope2:=Pekeris2(l,m,n,e,Z,L,M,N,3):
print(`We get that the coefficients`, A[l,m,n] , `are annihilated by the`):
print(`linear partial recurrence equation with polynomial coefficients`):
print(collect(ope2,[L,M,N])):

print(` where `, L,M,N, `are the shift operators in the discrete variables`):
print(l,m,n, `respectively. `):

lu:=0:

for i1 from ldegree(ope2,L) to degree(ope2,L) do
for i2 from ldegree(ope2,M) to degree(ope2,M) do
for i3 from ldegree(ope2,N) to degree(ope2,N) do
lu:=lu+
coeff(coeff(coeff(ope2,L,i1),M,i2),N,i3)*A[l+i1,m+i2,n+i3]:
od:
od:
od:
print(`Or more explicitly, in humanese`):
print(lu=0):

t1:=time():
print(`This ends the symbolic part, meticuously done by C.L.Pekeris`):
print(`It took`, t1-t0, `seconds `):
print(``):
print(`Now is time for the numeric part, done by WEIZAC and`):
print(`cleverly programmed in machine language by Yigal Accad`):
print(``):
print(`Taking account that we are talking about ortho states`):
print(`so that A[l,m,n]=-A[m,l,n], `):

print(`By truncating this infinite system to those`, l,m,n):
print(`satisfying `, l+m+n<=omega , `we get that the parameter`, e):
print(`where E=-e^2 ,for charges Z between 1, and `, MaxZ, `are as follows`):



for Z0 from 1 to MaxZ do
vu:=EnergyOrtho(omega,Z0,d,3,orekh):
print(`When the charge`,Z=Z0, `the appx. value is`,e=vu[1]):
print(`For omega=`,[seq(omega+i*d,i=-(orekh-1)..0)], `the e's for the truncated matrix are, resp.`):
print(vu[2]):
print(`The extrapolated values are`, vu[3]):
print(`the probable error is`, vu[4]):
od:

print(``):
print(`The whole thing took`, time()-t0, `seconds `):
end:
















#MatrixParaYigal(Z,e,om,memad): makes the set of equations for
#Pekeris's coeff.'s A(l,m,n) for 0<=l+m+n<=om
#in memad dimensions
#The output is the set of equations followed by the
#set of variables.For example, do:
#MatrixParaYigal(Z,e,1,3);
MatrixParaYigal:=proc(Z,e,om,memad) local l,m,n,l1,m1,n1,T1,ope,L,M,N,
lu,co,Mat1,ka1,ka2,l2,m2,n2,mu,a,b,i1,i2,i3,ku,Mat2,i:
option remember:
ope:=Pekeris2(l,m,n,e,Z,L,M,N,memad):
co:=0:

for n1 from 0 to om do
 for l1 from 0 to trunc((om-n1)/2) do
  for m1 from l1 to om-n1-l1 do
   co:=co+1:
   lu:=Klmn(l1,m1,n1):
    T1[lu]:=[l1,m1,n1]:
  od:
 od:
od:


for a from 1 to co do
 ka1:=T1[a]:  l1:=ka1[1]: m1:=ka1[2]: n1:=ka1[3]:
 for b from 1 to co do
   ka2:=T1[b]:  l2:=ka2[1]: m2:=ka2[2]: n2:=ka2[3]:
   i1:=l2-l1: i2:=m2-m1: i3:=n2-n1:
   mu:=subs({l=l1,m=m1,n=n1},coeff(coeff(coeff(ope,L,i1),M,i2),N,i3)):
    if l2<>m2 then
   i1:=m2-l1: i2:=l2-m1: i3:=n2-n1:
     mu:=mu+subs({l=l1,m=m1,n=n1},coeff(coeff(coeff(ope,L,i1),M,i2),N,i3)):
    fi:
   Mat1[a,b]:=mu:
 od:
od:

Mat1:=[seq([seq(Mat1[a,b],b=1..co)],a=1..co)]:

Mat1:

Mat2:=[]:

for i from 1 to nops(Mat1) do
 if Mat1[1][i]<>0 then
  ku:=normal(Mat1[i][1]/Mat1[1][i]):
 else
   ku:=1:
 fi:
 if not type(ku,numeric) then
   print(`Something is wrong`):
   RETURN(FAIL):
 fi:
 
   Mat2:=[op(Mat2),Mat1[i]/ku]:


od:

Mat2/(-2):


end: