function y2p = SecondDerivative(f,a,h)
% f is a handle to a function f(x) that has six continuous derivatives.
% y2p is an estimate of f"(a)

y2p = (-f(a-2*h) + 16*f(a-h) -30*f(a) + 16*f(a+h) - f(a+2*h))/(12*h^2);

% Error in y2p:
% f"(a) = y2p + (1/90)M6*h^4 M6 = f(6)(eta), a-2h <= eta <= a+2h

% Model overall error:
% err(h) = M6*h^4/90 + eps/h^2
% Taking derivatives:
% hopt = (45*eps/M6)^(1/6)

% Derivation of the recipe for y2p...

% f(a+h) = f(a) + hf'(a) + (h^2/2)f''(a) + (h^3/6)f'''(a) + (h^4/24)f""(a)
% + (h^5/120)f""'(a) + O(h^6)
% f(a-h) = f(a) - hf'(a) + (h^2/2)f''(a) - (h^3/6)f'''(a) + (h^4/24)f""(a)
% - (h^5/120)f""'(a) + O(h^6)
%
% (f(a-h)-2f(a)+f(a+h))/h^2 = f"(a) + (h^2/12)f""(a) + O(h^4)

% f"(a) = (f(a-h) -2f(a)+f(a+h))/h^2 - (h^2/12)f""(a) + O(h^4) (1)

% f"(a) = (f(a-2h)-2f(a)+f(a+2h))/(4h^2) - (h^2/3)f""(a) + O(h^4) (2)

% Multiply (1) by 4 and substract (2)...

% 3f"(a) = (4f(a-h)-8f(a)+4f(a+h))/h^2 - (f(a-2h)-2f(a)+f(a+2h))/(4h^2)
% = (4f(a-h)-8f(a)+4f(a+h) - f(a-2h)/4 + f(a)/2 -f(a+2h)/4 ) / h^2
% = (4f(a-h) + 4f(a+h) - f(a-2h)/4 - f(a+2h)/4 -(15/2)f(a))/h^2
% Thus,
% f"(a) = ((4/3)f(a-h) + (4/3)f(a+h)
% -(1/12)f(a-2h)-(1/12)f(a+2h)-(5/2)f(a))/h^2

% Partial Credit
% Correct coefficients: 3 points.
% Correct value of c: 1 point.
% Correct optimal value of h: 1 point.