function [D,x] = cheb(N)

% function [D,x] = cheb(N).
% Compute:
%   D = Chebyshev differentiation matrix
%   x = Chebyshev grid  
% given maximum polynomial degree N.
% From Trefethen, N, 2000: Spectral Methods in MATLAB, p. 54.

  if N==0, D = 0, x= 1; return, end
  x = cos(pi*(0:N)/N)';
  c = [2; ones(N-1,1); 2].*(-1).^(0:N)'; % c(j)*(-1)^j
  X = repmat(x,1,N+1);
  dX = X-X';                        % dX(i,j) = x(i) - x(j)
  D = (c*(1./c)')./(dX+eye(N+1));  % Off-diagonal entries 
                                    % D(i,j) = c(i)/c(j)*(-1)^(i-j)/(x(i)-x(j))
                                    % (eye() added to 
                                    % avoid division by zero on diagonal)
  D = D - diag(sum(D'));            % Round-off tolerant computation of 
                                    % diagonal elements based on zero row sums.
