% GRADAU_JACOBI Generalized Gauss-Radau quadrature formula
%    for Jacobi weight function
%
%    Given the parameters al and be of the Jacobi measure
%    and the fixed node -1 of multiplicity r, this routine
%    generates the generalized Gauss-Radau quadrature rule
%    involving r-1 consecutive derivative values at the
%    point -1 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 -1 (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, as it would occur in the routine gradau for
%    ab=r_jacobi.
%
function xw=gradau_jacobi(N,al,be,r)
%
% internal nodes and weights
%
apb=al+be+r;
ab1=r_jacobi(N,al,be+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)+1).^r);
xw(1:r,1)=-ones(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)+1).^(-rho));
end
p=1;
for n=1:N
  p=(1+(be+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_jacobi(ng,al,be);
xwg=gauss(ng,ab);
for i=1:r
  s=0;
  for ig=1:ng
    t=xwg(ig,1); w=xwg(ig,2);
    Pn=.5*(apb+2)*t+.5*(al-be-r); Pnm1=1;
    for n=2:N
      n2=2*n;
      Pnm2=Pnm1; Pnm1=Pn;
      Pn=((n2+apb-1)*((n2+apb)*(n2+apb-2)*t+al^2-(be+r)^2)*Pnm1 ...
        -2*(n+al-1)*(n+be+r-1)*(n2+apb)*Pnm2)/(n2*(n+apb)*(n2+apb-2));
    end
    s=s+w*((t+1)^(i-1)*Pn^2);
  end
  b(i)=s;
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;