format long
x = [0.6 0.7 0.8 0.9];
f = [cos(0.6) cos(0.7) cos(0.8) cos(0.9)];
y = f-x;
Q = zeros(4);
Q(:,1) = x';
for i = 1:3
  for j = 1:i
    Q(i+1,j+1) = ((0-y(i-j+1))*Q(i+1,j)-(0-y(i+1))*Q(i,j))/(y(i+1)-y(i-j+1));
    end
    end
Q(4,4)