function D = chebyshev_matrix ( n )

%*****************************************************************************80
%
%% chebyshev_matrix() computes the Chebyshev differentiation matrix.
%
%  Discussion:
%
%    Given a grid function v defined on the points of a Chebyshev grid,
%    the discrete derivative w is found by w = D * v.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 April 2020
%
%  Author:
%
%    Original MATLAB version by Lloyd Trefethen.
%    This version by John Burkardt.
%
%  Reference:
%
%    Lloyd Trefethen,
%    Spectral methods in MATLAB,
%    SIAM, 2000,
%    LC: QA377.T65
%    ISBN: 978-0-898714-65-4
%
%  Input:
%
%    integer n: the order of the grid.
%
%  Output:
%
%    real d(n+1,n+1): the differentiation matrix.
%
  if ( n == 0 )
    D = []; 
  else
%
%  Compute the Chebyshev points.
%
    x = cos ( pi * (0:n) / n )'; 

    c = [ 2.0; ones(n-1,1); 2.0 ] .* (-1.0) .^ (0:n)';
    X = repmat ( x, 1, n + 1 );
    dX = X - X';
%
%  Set the offdiagonal entries.
%       
    D  = ( c * ( 1.0 ./ c )' ) ./ ( dX + ( eye ( n + 1 ) ) );
%
%  Set diagonal entries.
%
    D  = D - diag ( sum ( D' ) );

  end

  return
end
