function ShowEveryThing()

close all
ShowPCHeur(.005,2000,4)
ShowPCHeur(.01,1000,4)
ShowPCHeur(.01,1000,2)
ShowPCHeur(.05,200,2)

function ShowPCHeur(h,n,d)
% Applies the d-th order Adams Predictor-Corrector Framework to
%
% y'(t) = f(t,y(t)) y(t0) = y0
%
% The steplength h is fixed and n steps are taken.
%
% Examines the heuristic that the discrepancy between the predicted and
% corrected values can be used to estimate the LTE.

% We use the rabbit-fox simulation.
f = @MyF;
t0 = 0;
rabbit0 = 300;
fox0 = 150;
y0 = [rabbit0;fox0];

% We need d f-evals to launch the d-th order AB method. Get them
% via the d-th order RK method...
[tvals,yvals,fvals] = ABStart(f,t0,y0,h,d);
tc = tvals(d);
yc = yvals(d,:)';

% For the LTE computation, we will use ode45 to compute the
% "true" yNext given that the current solution yc is correct.
options = odeset('Abstol',.0000000001);
% We will assemble the discrepancies between the predicted and
% corrected values in
Discrepancy = [];
% and the LTE in yCorr in
LTE= [];
% Right now we solutions f-evals at t = t0, t0+h,...,t0+(k-1)h
% Get solutions at t = t0+kh,...,t0+nh
for j=d:n
% Take a predictor-corrector step...
[tNext,yPred,fPred,yCorr,fCorr] = PCStep(f,tc,yc,fvals,h,d);
% Use ode45 to get the true yNext using yc as the initial condition
[t45,y45] = ode45(f,[tc tc+h],yc,options);
yNext = y45(length(t45),:)';
% Record the discrepancy and LTE estimate:
Discrepancy = [Discrepancy max(abs(yCorr-yPred))];
LTE = [LTE max(abs(yCorr-yNext))];
tvals = [tvals; tNext];
yvals = [yvals; yCorr'];
fvals = [fCorr fvals(:,1:d-1)];
tc = tNext;
yc = yCorr;
end
% Now lets look at the LTE estimation heuristic.
% Some known LTE facts about the AB and AM methods. For the
% k-th order AB and AM methods write
%
% y(tc+h) = yPred + c1*h^(d+1)*(some value of (d+1)-st derivative)
% y(tc+h) = yCorr + c2*h^(d+1)*(some value of (d+1)-st derivative)
%
% The c1 and c2 values are known. Here they are for d=1:5
c1 = [ 1/2 5/12 3/8 251/720 95/288];
c2 = [-1/2 -1/12 -1/24 -19/720 -3/160];
% If the (d+1)-st derivative is constant then
%
% yCorr -yPred = (c1-c2)*h^(d+1)*const
%
% and so
%
% y(tc+h) - yCorr = c2*h^(d+1)*const
% = ((c2/(c1-c2))*(c1-c2)*h^(d+1)*const
% = ((c2/(c1-c2))*(yCorr-yPred)
%
% So let's compare LTE and (c2(d)/(c1(d)-c2(d)))*Discrepancy..
factor = c2(d)/(c1(d)-c2(d));
figure
subplot(2,1,1)
plot(tvals(d:n),LTE,tvals(d:n),abs(factor*Discrepancy))
legend('LTE','Estimated LTE')
set(gca,'fontsize',14)
title(sprintf('d = %1d h = %6.3f',d,h),'fontsize',14)
xlabel('Time')
subplot(2,1,2)
plot(tvals,yvals(:,1),'r',tvals,yvals(:,2),'g')
legend('Rabbits','Foxes')
title('Rabbits and Foxes','fontsize',14)
xlabel('Time')
set(gca,'fontsize',14)

function yp = MyF(t,y)
% y1 = rabbits, y2 = foxes
alpha = .01;
yp = [2*y(1)-alpha*y(1)*y(2) ; -y(2)+alpha*y(1)*y(2)];