function e=h2plus_varzeta(zeta)
% to determine energy of 1sigma g state of H2+ using 1sig=(1sa+1sb)
% the necessary integrals are given on page 6 of Chapter 3
   global r
   rho=zeta*r;
   rho2=rho.*rho;
   s=(1+rho+rho2/3);
   ex=exp(-rho);
   s=s.*ex;
   ke=-0.5*zeta*zeta*(-1-rho+rho2/3).*ex;
   vab=-zeta*(1+rho).*ex;
   vaa=-zeta*(1-(1+rho).*ex.*ex)./rho;
   e=(zeta^2-2*zeta+2*ke+4*vab+2*vaa)./(2*(1+s));