clear
clf

h = 1;
np = 16;
x = 1.;

hh = zeros(np, 1);
der2a = zeros(np, 1);
der2b = zeros(np, 1);

f = @(x) sin(x);

for i = 1:np
    hh(i) = h;
    der2a(i) = (f(x - h) - 2*f(x) + f(x + h))/(h^2);
    der2b(i) = (-f(x - 2*h) + 16*f(x - h) - 30*f(x) +  16*f(x + h) - f(x + 2*h))/(12*h^2);
    h = h/sqrt(2);
end

loglog(hh, abs(der2b + sin(x)), 'bo-')
hold on
grid on
loglog(hh, hh.^4/1000, 'b--')

loglog(hh, abs(der2a + sin(x)), 'ro-')
loglog(hh, hh.^2, 'r--')

xlabel('h');
ylabel('\delta(h)');

legend('4th order','h^4','2nd order','h^2','location','best')
title('Accuracy of finite difference derivatives')