x = 1;

f = @(y) sin(y);

h = 1;
np = 15;

sq = sqrt(2.);

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

for i = 0:np
    hh(i+1) = h;
    der2(i+1) = (-f(x-2*h) + 16*f(x-h) -30*f(x) + 16*f(x+h) - f(x+2*h))/(12*h^2);
    h = h/sq;
end

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