format long

f = [0.3, 0.4, 0.5, 0.6];
x = [0.5488, 0.4493, 0.3679, 0.3012];
x = f - x;
xt = 0;
d = length(x);
Q = zeros(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)