% D=lafgsrddiff(n,x) returns the first-order differentiation matrix (of the Laguerre function approach) of size
% n by n, at the Laguerre-Gauss-Radau points x, which can be computed by 
% x=lagsrd(n). Note: x(1)=0.
% See Page 259 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: lapoly() 
% Last modified on December 22, 2011

function D=lafgsrddiff(n,x)
if n==0, D=[]; return; end;
nn=n-1;
y=lafun(nn,x);  % Compute L_{n-1}(x_j) 
nx=size(x); xx=x; ny=size(y); 
if nx(1)>nx(2), xx=x'; end;  %%  xx is a row vector 
if ny(2)>nx(1), y=y'; end;   %% y is a column vector
  xx=xx(2:end); y=y(2:end);         % take x_k: k=2,..., n; note: x_1=0.
  
  %Row 1:
  d1=[-nn/2,-1./(xx.*y')]; 
  %Row 2 to n
    %% Column 1
    d2=y./xx'; 
    %% Column 2 to n 
    D=(xx'./y)*y'-(1./y)*(xx.*y');  %% compute dy(x_j) (x_k-x_j)/dy(x_k);     
                                   % 1/d_{kj} for k not= j (see (7.49)) 
    D=D+eye(nn);                    % add the identity matrix so that 1./D can be operated                                     
    D=1./D; D=D-eye(nn);   D=D+eye(nn)/2 ;  %update the diagonal entries  
    
  % Assemble the whole matrix  for Lagueree polynomial approach
   D=[d1;[d2,D]];  
   D=D-diag(0.5*ones(size(x))); %% by (7.52);
   
  return;