% Problem specification.
f = @(x) sin(x);
x0 = 0.5;

% Exact solution (for comparision).
format longe
Q = cos(x0)

% Forward difference rule.
F = @(h) (f(x0 + h) - f(x0)) / h;

h = 0.1;  % Initial stepsize
M = 12;   % Number of iterations

% The extrapolation process, based on the error expansion
% F(h) = Q + C1*h^1 + C2*h^2 + C3*h^3 + ...
% and with "harmonic" reduction of the stepsizes: h / j, with j = 1, 2, ...

T = zeros(M);  
for j = 1:M
    T(j, 1) = F(h / j);
    for k=1:j-1
        T(j, k+1) = T(j, k) + (T(j, k) - T(j-1, k)) / ((j / (j - k)) - 1);
    end  
end

% Display extrapolation table and errors.
T
format shorte
Error = tril(Q-T)