function [x, DM] = herdif(N, M, b);
% The function [x, DM] = herdif(N, M, b) computes the
% differentiation matrices D1, D2, ..., DM on Hermite points.
%
% Input:
% N: Number of points, i.e., order of differentiation matrices (integer).
% M: Number of derivatives required (integer).
% b: Scaling parameter (real, positive).
%
% Note: 0 < M < N-1.
%
% Output:
% x: Vector of nodes (zeros of Hermite polynomial of degree N,
% scaled by the parameter b.)
% DM: DM(1:N,1:N,l) contains l-th derivative matrix, l=1..M.
% J.A.C Weideman, S.C. Reddy 1998.
x = herroots(N); % Compute Hermite roots.
alpha = exp(-x.^2/2); % Compute weights.
beta(1,:) = ones(size(x')); % Set up beta matrix s.t. beta(l,j) =
beta(2,:) = -x'; % (l-th derivative of alpha(x))/alpha(x),
% evaluated at x = x(j).
for ell = 3:M+1
beta(ell,:) = -x'.*beta(ell-1,:)-(ell-2)*beta(ell-2,:);
end
beta(1,:) = []; % Remove initializing row from beta
DM = poldif(x, alpha, beta); % Compute differentiation matrix (b=1).
x = x/b; % Scale nodes by the factor b.
for ell = 1:M % Adjust for b not equal to 1.
DM(:,:,ell) = (b^ell)*DM(:,:,ell);
end
