%Example of using Richardson Extrapolation when the error is
% of the form:  F(h)+k_1 h^2 +k_2 h^4 + k_3 h^6 + ...
%
%  In this case, use F(h)=f'(2)=(f(2+h)-f(2-h))/2h
%  using f(x)=xexp(x)  (Actual value:  3exp(2))


f=inline('x.*exp(x)');
h=0.2;
x=2;

for j=1:4
  H=h/2^(j-1);
  F(j)=(f(x+H)-f(x-H))/(2*H);
end
n=length(F);
Q=zeros(n,n);
Q(:,1)=F(:);

for col=2:n
  for row=col:n
    Q(row,col)=(4^(col-1)*Q(row,col-1)-Q(row-1,col-1))/(4^(col-1)-1);
  end
end
%The solution is 3e^2, so subtract that from everything to see the errors:
Error=diag(abs(Q-3*exp(2)));


% OUTPUT should be:
%Error =
%
%   0.246992360237467
%   0.000172675392047
%   0.000000013206471
%   0.000000000000227