function P2()
% Checks out PeriodicSimpson
a = 1;
P = 10; % The period
M4 = 128*(pi/P)^4; % A bound for the 4th derivative
tol = .000000001; % Absolute error tolerance
clc
fprintf(' k tau |I = PerSimp| f-evaluations\n')
fprintf('-------------------------------------------------------\n')
for k=[0 10 1000]
% The number of full periods
for j=1:10
tau = j/10; % The fractional period
% Integrate MyF from a to b where
b = a + (k+tau)*P;
[Q,fEvals] = PeriodicSimpson(@(x) MyF(x,P),a,b,P,M4,tol);
I = (b-a)/2 + (P/(8*pi))*(sin(4*pi*b/P) - sin(4*pi*a/P));
err = abs(Q - I);
fprintf(' %4d %5.1f %10.3e %6d \n',k,tau,err,fEvals)
end
end

function y = MyF(x,P)
y = cos(2*pi*x/P).^2;
% The integral of this from a to b is
% I = (b-a)/2 + (P/8*pi)*(sin(4*pi*b/P) - sin(4*pi*a/P))
% The 4th derivative bound is 128*(pi/P)^4