function d = div_dif(x,f)

% function d = div_dif(x,f)
% divided differences d_k = Delta(x_1,...,x_k) f 
% based on Taylor coefficients f_j, f'_j, f''_j/2, ... 
% according to sequence of points x_i

n = length(x);
for k=1:n-1;
   dx = x(k+1:n)-x(1:n-k);
   ind = find(dx>0);
   g = [f(k); f(k+ind)];
   f(k+ind) = diff(g)./dx(ind);
end;
d = f;

