clear
clf

h = 1;
np = 16;
x = 1.;
sq = sqrt(2);

hh = zeros(np, 1);
errder1a = zeros(np, 1);
errder1b = zeros(np, 1);

f = @(x) sin(x);

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

loglog(hh, abs(errder1b), 'b.-')
xlabel('h')
ylabel('\delta(h)')
title('Accuracy of finite difference derivatives')
hold on
loglog(hh, hh.^4/1000, 'b--')

loglog(hh, abs(errder1a), 'r.-')
loglog(hh, hh.^2, 'r--')

legend('4th order','h^4','2nd order','h^2','location','best')
grid on