function F=divdiff(x,y)
n=length(x)-1;
for i=1:n+1
    F(i,1)=y(i);
end
for i=2:n+1
    for j=2:i
        F(i,j)=(F(i,j-1)-F(i-1,j-1))/(x(i)-x(i-j+1));
    end
end

% d=divdiff([1.0 1.3 1.6 1.9 2.2],[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623])
% save 'divdiff.txt' d -ascii