function [x, DM] = sincdif(N, M, h)%  The function [x, DM] = sincdif(N, M, b) computes sinc the%  differentiation matrices D1, D2, ..., DM on equidistant points.%%  Input:%  N:    Number of points, i.e., order of differentiation matrix.%  M:    Number of derivatives required (integer).%  h:    Step-size (real, positive).%%  Note:  0 < M < N-1.%%  Output:%  x:    Vector of nodes.%  DM:   DM(1:N,1:N,l) contains l-th derivative matrix, l=1..M.%  J.A.C. Weideman, S.C. Reddy 1998.k = [1:N-1]';t = k*pi;x = h*[-(N-1)/2:(N-1)/2]';sigma = zeros(size(k));for ell = 1:M;    sigma = (-ell*sigma + imag(exp(i*t)*i^ell))./t;      col = (pi/h)^ell*[imag(i^(ell+1))/(ell+1); sigma];      row = (-1)^ell*col; row(1) = col(1);DM(:,:,ell) = toeplitz(col,row);end