function P1(delta,h)
% Use AM2 methods to solve
%
% y'(t) = y(t)^2 - y(t)^3, y(0) = delta
%
% across [0,nh] where
n = ceil((2/delta)/h);

% Where we approximate y(t):
tFinal = n*h;
tVals = linspace(0,tFinal,n+1)';

% Initial time and initial condition...
t0 = 0;
y0 = delta;

% Use ode23s and spline to get the "exact" solution...
opts = odeset('RelTol',.000001);
[tOde23s,yOde23s] = ode23s(@flame,[t0 tFinal],y0,opts);
yExact = spline(tOde23s,yOde23s,tVals);

% Use the fixed point-based method...
close all
global fCount
fCount=0;
yOut = AM2FixPt(@flame,t0,y0,n,h);
plot(tVals,yOut)
title(sprintf('Fixed Point-Based Method: h = %6.4f fEvals = %1d Error = %10.3e',h,fCount,max(abs(yOut-yExact))))

% Use the Newton-based method...
fCount=0;
yOut = AM2Newton(@flame,t0,y0,n,h);
figure
plot(tVals,yOut)
title(sprintf('Newton-Based Method: h = %6.4f fEvals = %1d Error = %10.3e',h,fCount,max(abs(yOut-yExact))))
shg

function [yp,dfdy] = flame(t,y)
global fCount
yp = y^2-y^3;
if nargout==2
dfdy = 2*y - 3*y^2;
end
fCount = fCount+1;