% QUADSQEINSTEIN A quadrature routine used in Example 2.42.
%
function xw=quadsqeinstein(N,i)
global r
c1=10*sqrt(10); c2=exp(-10*r)/r;
if i==1
  ab=r_jacobi01(N,0,1/2); xw=gauss(N,ab);
  xw(:,1)=10*xw(:,1);
  for n=1:N
    t=xw(n,1);
    if abs(t)<1
      term=1; s=0; s1=1; k=1;
      while s1~=s
        s=s1;
        k=k+1;
        term=-t*term/k; s1=s+term;
      end
      xw(n,2)=c1*xw(n,2)*exp(-r*t)/(s^r);
    else
      xw(n,2)=c1*xw(n,2)*(t/(exp(t)-1))^r;
    end
  end
else
  ab=r_laguerre(N); xw=gauss(N,ab);
  t=10+xw(:,1)/r;
  xw(:,1)=t;
  xw(:,2)=c2*xw(:,2).*sqrt(t).*(t./(1-exp(-t))).^r;
end
