% D=lagsdiff(n,x) returns the first-order differentiation matrix of size
% n by n, associated with the Laguerre-Gauss points x, which may be computed by 
% x=lags(n). 
% See Page 251 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=lagsdiff(n,x)
if n==0, D=[]; return; end;
xx=x;y=lapoly(n-1,xx); nx=size(x); 
if nx(2)>nx(1), y=y'; xx=x'; end;  %% xx, y are column-n vectors
  D=(xx.^2./y)*(y./xx)'-(xx./y)*y';  %% 
  D=D+eye(n);                    % add the identity matrix so that 1./D can be operated                                     
  D=1./D; 
  D=D-eye(n); D=D+diag(1/2-1./(2*xx'));  % update the diagonal entries  
      
  return;