  function Aug26()
% f(a+h) = f(a) + h*f'(a) + (h^2/2)f''(a) + (h^3/6)f'''(a) + (h^4/24)f''''(a) + O(h^5)
% f(a-h) = f(a) - h*f'(a) + (h^2/2)f''(a) - (h^3/6)f'''(a) + (h^4/24)f''''(a) + O(h^5)

% So
%        f(a+h)-f(a-h)   = 2hf'(a) + (h^3/3)f'''(a) + O(h^5)    (1)

% The same formula with h replaced by 2h...

%        f(a+2h)-f(a-2h) = 4hf'(a) + (8h^3)f'''(a) + O(h^5)     (2)
%
% Multiply (1) by 8 and subtract (2)...
% 
%     8f(a+h)-8f(a-h)-f(a+2h)+f(a-2h) = 12hf'(a) + O(h^5)
%
% so
%    f'(a) = (8f(a+h)-8f(a-h)-f(a+2h)+f(a-2h))/12h  + O(h^4)

clc
disp('    h         error ')
disp('----------------------------')
for k=1:10
    h = 1/2^k;
    yp = Derivative4(@sin,1,h);
    err = abs(yp-cos(1));
    fprintf('%10.3e  %10.3e \n',h,err)
end

function yp = Derivative4(f,a,h)
yp = (8*f(a+h)-8*f(a-h)-f(a+2*h)+f(a-2*h))/(12*h);
