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
