function fn=divdif(xn,yn)
   fn=yn;
   n=length(xn);
   for k=1:(n-1)
     for i=n:-1:(1+k)
       fn(i)=(fn(i-1)-fn(i))/(xn(i-k)-xn(i));
     end
   end