function P2()

clc
% Apply SecondDerivative to compute the 2nd derivative of x^6 at x = a;
% The true answer is 30.
% The sixth derivative of this function is the constant 6! = 720.
a = 1;
for h = logspace(-1,-2,20)
y2p = SecondDerivative(@(x) x^6,a,h);
% 30*a^4 = y2p + const*(720)*h^4
const = (30*a^4 - y2p)/(720*h^4);
fprintf('%10.3e %10.6f\n',h,const)
end

% From which we deduce that the const = .011111... = 1/90

% Let's try SecondDerivative out on two examples
close all
% Example 1. f(x) = sin(x), estimate f"(1)
a = 1;
N = 100;
h = logspace(-5,-1,N);
for k=1:N
y2p(k) = SecondDerivative(@sin,a,h(k));
end
loglog(h,abs(y2p+sin(a))+eps)
hold on
M6 = 1;
hopt = (45*eps/M6)^(1/6);
plot([hopt hopt],ylim,'r')
hold off
title(sprintf('f(x) = sin(x), a = %5.3f, M6 = 1, hopt = %6.4f',a,hopt),'fontsize',12)
xlabel('h','fontsize',12)
ylabel('| f"(a) - SecondDervative(f,a,h) |','fontsize',12)

% Example 2. f(x) = sin(10x), estimate f"(1)
figure
a = 1;
N = 100;
h = logspace(-5,-1,N);
for k=1:N
y2p(k) = SecondDerivative(@(x) sin(10*x),a,h(k));
end
loglog(h,abs(y2p+100*sin(10*a))+eps)
hold on
M6 = 1000000;
hopt = (45*eps/M6)^(1/6);
plot([hopt hopt],ylim,'r')
hold off
title(sprintf('f(x) = sin(10x), a = %5.3f, M6 = 1000000, hopt = %6.4f',a,hopt),'fontsize',12)
xlabel('h','fontsize',12)
ylabel('| f"(a) - SecondDervative(f,a,h) |','fontsize',12)
shg