function [fx,dfx] = hornerNewton(a,x,x0)
% evaluate Newton polynomial and its derivative
% using Horner's scheme
% p(x0) = a(1) + a(2)*(x0-x(1)) + a(3)*(x0-x(1))*(x0-x(2)) + ...
%       = a(1) + ( (x0-x(1)) * ( a(2) + (x0-x(2)) * ( a(3) .... 
fx = a(end);  dfx = 0;
for k=length(a)-1:-1:1
  dfx = fx + (x0-x(k)).*dfx;
  fx =  a(k) + (x0-x(k)).*fx; 
end

