function y=ddfun(X,Y,x)
F=divdiff(X,Y);
n=length(X);
y=0;
for i=1:n
    p=F(i,i);
    for j=1:i-1
        p=p.*(x-X(j));
    end
    y=y+p;
end
       
% ddfun([1.0 1.3 1.6 1.9 2.2],[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623],1.5)