x = 1;

f = @(y) sin(y);

h = 1;
np = 15;

sq = sqrt(2.);

hh = zeros(np+1, 1);
der1 = zeros(np+1, 1);

for i = 0:np
    hh(i+1) = h;
    der1(i+1) = (f(x-2*h) - 8*f(x-h)  + 8*f(x+h) - f(x+2*h))/(12*h);
    h = h/sq;
end

loglog(hh, abs(der1 - cos(x)), 'ro-');
hold on
loglog(hh, hh.^4, 'b--');
grid on
% Add labels, title, etc.
