function [x, DM] = chebdif(N, M)%  The function DM =  chebdif(N,M) computes the differentiation %  matrices D1, D2, ..., DM on Chebyshev nodes. % %  Input:%  N:        Size of differentiation matrix.        %  M:        Number of derivatives required (integer).%  Note:     0 < M <= N-1.%%  Output:%  DM:       DM(1:N,1:N,ell) contains ell-th derivative matrix, ell=1..M.%%  The code implements two strategies for enhanced %  accuracy suggested by W. Don and S. Solomonoff in %  SIAM J. Sci. Comp. Vol. 6, pp. 1253--1268 (1994).%  The two strategies are (a) the use of trigonometric %  identities to avoid the computation of differences %  x(k)-x(j) and (b) the use of the "flipping trick"%  which is necessary since sin t can be computed to high%  relative precision when t is small whereas sin (pi-t) cannot.    %  J.A.C. Weideman, S.C. Reddy 1998.     I = eye(N);                          % Identity matrix.          L = logical(I);                      % Logical identity matrix.    n1 = floor(N/2); n2  = ceil(N/2);     % Indices used for flipping trick.     k = [0:N-1]';                        % Compute theta vector.    th = k*pi/(N-1);     x = sin(pi*[N-1:-2:1-N]'/(2*(N-1))); % Compute Chebyshev points.     T = repmat(th/2,1,N);                    DX = 2*sin(T'+T).*sin(T'-T);          % Trigonometric identity.     DX = [DX(1:n1,:); -flipud(fliplr(DX(1:n2,:)))];   % Flipping trick.  DX(L) = ones(N,1);                       % Put 1's on the main diagonal of DX.     C = toeplitz((-1).^k);               % C is the matrix with C(1,:) = C(1,:)*2; C(N,:) = C(N,:)*2;     % entries c(k)/c(j)C(:,1) = C(:,1)/2; C(:,N) = C(:,N)/2;     Z = 1./DX;                           % Z contains entries 1/(x(k)-x(j))    Z(L) = zeros(N,1);                      % with zeros on the diagonal.     D = eye(N);                          % D contains diff. matrices.                                          for ell = 1:M          D = ell*Z.*(C.*repmat(diag(D),1,N) - D); % Off-diagonals       D(L) = -sum(D');                            % Correct main diagonal of DDM(:,:,ell) = D;                                   % Store current D in DMend