function yOut = AM2FixPt(f,t0,y0,n,h)
% f is a handle to a function of the form yp = f(t,y) where t and y are scalars.
% yOut is a column (n+1)-vector with the property that
% yOut(k) is an approximate solution to the IVP
%
% y'(t) = f(t,y) y(t0) = y0
%
% at t = t0+(k-1)h, k=1:n+1. The 2nd order Adams-Moulton method (AM2)
% with step length h is used to generate the approximate solutions. In
% particular, z = y_{k+1} approximately satisfies
%
% z = y_k + (h/2)*(f(t_{k+1},z) + f(t_{k},y_{k})
%
% where t_{k}=t0+kh and y_{k} \approx y(t_[k}).
% Fixed point iteration is used to do this.

yOut = zeros(n+1,1);
% Since the local truncation error is O(h^3), there is no point
% point in trying to satisfy the AM2 equation to greater accuracy.
tol = h^3;
tc = t0;
yOut(1) = y0;
fc = f(tc,yOut(1));
fcount = 0;
for k=1:n
[yOut(k+1),fc] = AM2FixPtStep(f,tc,yOut(k),fc,h,tol);
tc = tc+h;
end

function [yNext,fNext] = AM2FixPtStep(f,tc,yc,fc,h,tol)
% Integrates y' = f(t,y) from tc to tNext = tc+h assuming y(tc) = yc
% and fc = f(tc,yc).
% Uses the second-order Adams-Moulton method.
% yNext approximately satisfies the equation
%
% yNext = yc + (h/2)(f(tnext,yNext) + fc);
%
% The function evaluation f(tnext,yNext) is returned in fNext.

% Find a fixed point of G(z) = (h/2)*f(tNext,z) + const where
tNext = tc+h;
const = yc + (h/2)*fc;
j=1;
% zc and zNext house the (j-1)-st and jth iterate. Our initial guess
% is the approximate current solution yc since it is "free".
zc = yc;
zNext = (h/2)*fc + const;
fNext = f(tNext,zNext);
while abs(zc-zNext)> tol && j<100
zc = zNext;
zNext = (h/2)*fNext + const;
fNext = f(tNext,zNext);
j = j+1;
end
yNext = zNext;