% GRADAU_LAGUERRE Generalized Gauss-Radau quadrature formula
%    for generalized Laguerre weight function
%
%    Given the parameter al of the generalized Laguerre weight
%    function and the fixed node 0 of multiplicity r, this routine
%    generates the generalized Gauss-Radau quadrature rule
%    involving r-1 consecutive derivative values at the
%    point 0 and N interior points. The N interior nodes
%    and weights are stored in the last N rows of the
%    (N+r)x2 output array xw, and the fixed node 0 (repeated
%    r times) and the r weights associated with the r
%    derivative values in the first r rows of xw.
%
%    The routine is designed to avoid breakdown at very large
%    values of N, which would occur in the routine gradau for
%    ab=r_laguerre. 
%
function xw=gradau_laguerre(N,al,r)
%
% internal nodes and weights
%
ab1=r_laguerre(N,al+r);
xw0=gauss(N,ab1);
xw(r+1:N+r,:)=xw0;
xw(r+1:N+r,2)=xw(r+1:N+r,2)./(xw(r+1:N+r,1).^r);
xw(1:r,1)=zeros(r,1);
%
% boundary weights
%
rs=zeros(r-1,1); 
A=zeros(r); b=zeros(r,1); x=zeros(r,1);
for rho=1:r-1
  rs(rho)=sum(xw0(:,1).^(-rho));
end
p=1;
for n=1:N
  p=(1+(al+r)/n)*p;
end
for i=r:-1:1
  A(i,i)=factorial(i-1)*p^2;
  if i<r
    for j=i+1:r
      A(i,j)=-2*sum(A(i+1:j,j).*rs(1:j-i))/(j-i);
    end
  end
end
condA=cond(A);
ng=N+floor((r+1)/2);
ab=r_laguerre(ng,al);
xwg=gauss(ng,ab);
for i=1:r
  s0=0; ratn=1;
  for ig=1:ng
    t=xwg(ig,1); w=xwg(ig,2);
    c=w*t^(i-1);
    Pn=c*(al+r+1-t); Pnm1=c;
    for n=2:N
      Pnm2=Pnm1; Pnm1=Pn;
      Pn=((2*n+al+r-1-t)*Pnm1-(n+al+r-1)*Pnm2)/n;
    end
    if c~=0
      s1=s0+Pn^2/c; 
      s0=s1;
    end
  end
  b(i)=s1;
end
x(r)=b(r)/A(r,r);
for i=r-1:-1:1
  x(i)=(b(i)-sum(A(i,i+1:r).*(x(i+1:r))'))/A(i,i);
end
xw(1:r,2)=x;