format long

f = [0.6, 0.7, 0.8, 0.9];
x = [cos(0.6)-0.6, cos(0.7)-0.7, cos(0.8)-0.8, cos(0.9)-0.9];

xt = 0;
d = length(x);
n = d + 1;
Q = ones(d);

for i = 1:d
  for j = 1:i
    if j == 1
      Q(i, j) = f(i);
    else
      Q(i, j) = (xt - x(i - j + 1))*Q(i, j - 1) - (xt - x(i))*Q(i - 1, j - 1);
      Q(i, j) = Q(i, j) / (x(i) - x(i - j + 1));
    end
      
  end
end
Q(d, d)