x = -1:0.01:1;
f = exp(x)-1;
y = f-x;
n = size(x,2);
Q = zeros(n);
Q(:,1) = x';
for i = 1:n-1
  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(n,n)