% QUADEINSTEIN A quadrature routine used in Exanple 2.41.
%
function xw=quadeinstein(N,i)
global r
ab=r_laguerre(N); xw=gauss(N,ab);
xw(:,1)=xw(:,1)/r;
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)=xw(n,2)/(r*s^r);
  else
    xw(n,2)=xw(n,2)*((t/(1-exp(-t)))^r)/r;
  end
end