function DM = poldif(x, malpha, B)%  The function DM =  poldif(x, maplha, B) computes the%  differentiation matrices D1, D2, ..., DM on arbitrary nodes.%%  The function is called with either two or three input arguments.%  If two input arguments are supplied, the weight function is assumed %  to be constant.   If three arguments are supplied, the weights should %  be defined as the second and third arguments.%%  Input (constant weight):%%  x:        Vector of N distinct nodes.%  malpha:   M, the number of derivatives required (integer).%  B:        Omitted.%%  Note:     0 < M < N-1.%%  Input (non-constant weight):%%  x:        Vector of N distinct nodes.%  malpha:   Vector of weight values alpha(x), evaluated at x = x(k).%  B:        Matrix of size M x N,  where M is the highest %            derivative required.  It should contain the quantities %            B(ell,j) = beta(ell,j) = (ell-th derivative%            of alpha(x))/alpha(x),   evaluated at x = x(j).%%  Output:%  DM:       DM(1:N,1:N,ell) contains ell-th derivative matrix, ell=1..M.%  J.A.C. Weideman, S.C. Reddy 1998       N = length(x);                             x = x(:);                     % Make sure x is a column vector if nargin == 2                       % Check if constant weight function       M = malpha;                   % is to be assumed.   alpha = ones(N,1);                     B = zeros(M,N);elseif nargin == 3   alpha = malpha(:);                % Make sure alpha is a column vector       M = length(B(:,1));           % First dimension of B is the number end                                  % of derivative matrices to be computed         I = eye(N);                  % Identity matrix.        L = logical(I);              % Logical identity matrix.       XX = x(:,ones(1,N));       DX = XX-XX';                  % DX contains entries x(k)-x(j).    DX(L) = ones(N,1);               % Put 1's one the main diagonal.        c = alpha.*prod(DX,2);       % Quantities c(j).        C = c(:,ones(1,N));         C = C./C';                   % Matrix with entries c(k)/c(j).           Z = 1./DX;                   % Z contains entries 1/(x(k)-x(j))     Z(L) = zeros(N,1);              % with zeros on the diagonal.        X = Z';                      % X is same as Z', but with      X(L) = [];                      % diagonal entries removed.        X = reshape(X,N-1,N);        Y = ones(N-1,N);             % Initialize Y and D matrices.        D = eye(N);                  % Y is matrix of cumulative sums,                                     % D differentiation matrices.for ell = 1:M        Y   = cumsum([B(ell,:); ell*Y(1:N-1,:).*X]); % Diagonals        D   = ell*Z.*(C.*repmat(diag(D),1,N) - D);   % Off-diagonals     D(L)   = Y(N,:);                                % Correct the diagonalDM(:,:,ell) = D;                                     % Store the current Dend