% D=legsrddiff(n,x) returns the first-order differentiation matrix of size
% n by n, associated with the Legendre-Gauss-Radau points x, which may be computed by 
% x=legsrd(n). Note: x(1)=-1.
% See Page 110 of the book: J. Shen, T. Tang and L. Wang, Spectral Methods:
%  Algorithms, Analysis and Applications, Springer Series in Compuational
%  Mathematics, 41, Springer, 2011. 
%  Use the function: lepoly() 
% Last modified on August 31, 2011

function D=legsrddiff(n,x)
if n==0, D=[]; return; end;
xx=x; nx=size(x);
[dy1,y1]=lepoly(n-1,xx); [dy2,y2]=lepoly(n,xx); dy=dy1+dy2; y=y1+y2; % Compute L_{n-1}+L_n and its 1st-order derivative
if nx(2)>nx(1), y=y'; dy=dy'; xx=x'; end;  %% y is a column vector of L_{n-1}(x_k)
  D=(xx./dy)*dy'-(1./dy)*(xx.*dy)';  %% compute dy(x_j) (x_k-x_j)/dy(x_k);     
                                 % 1/d_{kj} for k not= j (see (3.204)) 
  D=D+eye(n);                    % add the identity matrix so that 1./D can be operated                                     
  D=1./D; 
  D=D-eye(n); xx=xx(2:end);
  D=D+diag([-(n+1)*(n-1)/4; xx./(1-xx.^2)+ n*y1(2:end)./((1-xx.^2).*dy(2:end))]);  %update the diagonal entries  
  return;