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 diagonal
DM(:,:,ell) = D;                                     % Store the current D
end