function [A,df] = diffnew(X,Y)
% [A,df] = diffnew(X,Y)
% Numerical approximation for f'(x).
% The method is differentiation of the Newton polynomial Pn(x).
% An approximation for f'(X(1)) is computed.
% X is the list of abscissas, input.
% Y is the list of abscissas, input.
% A is the coefficient list for Pn(x), output.
% df is the approximate derivative, output.
A = Y;
n = length(X);
for j=2:n,
for k=n:-1:j,
A(k) = (A(k)-A(k-1))/(X(k)-X(k-j+1));
end
end
x0 = X(1);
df = A(2);
prod = 1;
n1 = length(A)-1;
for k=2:n1,
prod = prod*(x0 - X(k));
df = df + prod*A(k+1);
end