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

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

% Basic algorithm
F = @(h) (f(x0+h)-f(x0-h))/(2*h)

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

% The extrapolation process, based on the error expansion
% F(h) = Q + C1*h^2 + C2*h^4 + C3*h^6 + .....
% and by halving the stepsize h. 

T = [];  
for j = 1:M
    T(j,1) = F(h/2^(j-1));
    for k=1:j-1
        T(j,k+1) = (4^k*T(j,k)-T(j-1,k))/(4^k-1);
    end  
end

format long
T
format short e
Error = tril(Q-T)


