  function yOut = AM2Newton(f,t0,y0,n,h)
% f is a handle to a function of the form [yp,dfdt] = f(t,y) where t and y are scalars.
% yp is the rate of change of f(t,y(t)) with respect to t and dfdy is the
% rate of change of f with respect to y.
% 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}).
% A Newton 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;
  [fc,dfdy] = f(tc,y0);
  yOut(1) = y0;
  
  for k=1:n
     [yOut(k+1),fc,dfdy] = AM2NewtonStep(f,yOut(k),tc,fc,dfdy,h,tol);
     tc = tc+h;
  end
  
  


  function [yNew,fNew,dfdyNew] = AM2NewtonStep(f,yc,tc,fc,dfdy,h,tol) 
% Integrates y' = f(t,y) from tc to tc+h  assuming y(tc) = yc
% and fc = f(tc,yc); Uses the second-order Adams-Moulton method.
% ynew satisfies the equation
%
%           ynew = yc + (h/2)(f(tnew,ynew) + fc);
%
% where tnew = tc+h. The function evaluation f(tnew,ynew) is returned in
% fnew.

% Apply Newton to find a zero of
%
%                  F(z) = z - (yc + (h/2)*(f(tnew,z)+fc)
% i.e.,
%                  F(z) = z - (h/2)f(tNew,z) + const
% where

const = -yc - (h/2)*fc;

% Note  that F'(z) = 1 - (h/2)*dfdy(tnew,z)
% Our initial guess
% is the approximate current solution yc since it is "free".
z = yc;
F = z - (yc + h*fc);
tNew = tc+h;
j = 0;
while j==0 || (abs(F) >= tol && j<=10)
    Fp = 1 - (h/2)*dfdy;
    z = z - F/Fp;
    [fNew,dfdy] = f(tNew,z);
    F = z - (h/2)*fNew + const;
    j = j+1;
end
yNew = z;
dfdyNew = dfdy;