f=inline('x.*exp(x)');
h=0.1;
x=1;

for j=1:9
  H=h/2^(j-1);
  F(j)=(f(x+H)-f(x))/H;
end
n=length(F);
Q=zeros(n,n);
Q(:,1)=F(:);
for col=2:n
  for row=col:n
    Q(row,col)=(2^(col-1)*Q(row,col-1)-Q(row-1,col-1))/(2^(col-1)-1);
  end
end
%The solution is 2e, so subtract that from everything to see the errors:
Error=diag(abs(Q-2*exp(1)));

plot(log10(Error));
title('Log (base 10) of the error versus extrapolation number')