% the program takes the initial data y_0 at the initial time x_0 and 
% computes up to time x_final using N timesteps of size dx.  It creates
% two vectors, x and y, which you can then plot against each other
% plot(x,y)
%
% To call the routine, use the command:
%
% [x,y] = adams_moulton_2(y_0,x_0,x_final,N,@f)
%
function [x,y] = adams_moulton_2(y_0,x_0,x_final,N,f)

% step-size:
dx = (x_final-x_0)/N;

% for dy/dx = f(x,y(x))
%
% y(x+dx) ~ y(x) + dx/2*(f(x,y(x)) + f(x+dx,y(x+dx)))
%
% AKA Crank-Nicolson.  This is implicitly determining y(x+dx)
% so we will do an approximation of y(x+dx) which will then
% go into the right-hand side...

% pux initial data into vector
x(1) = x_0;
y(:,1) = y_0;
% I will take one Euler step.  This has local truncation error O(h^2).
F_now = feval(f,x(1),y(:,1));
x(2) = x_0 + dx;
y(:,2) = y(:,1) + dx*F_now;
for ii=2:N
  F_old = F_now;
  F_now = feval(f,x(ii),y(:,ii));
  % I make a first approximation of y(x+dx) using the second-order
  % Adams-Bashforth
  y_approx = y(:,ii) + dx/2*(3*F_now - F_old);
  % I make a second approximation of y(x+dx) as if finding the root
  % of z = y(x) + dx/2*(f(x,y) + f(x+dx,z)) using an iterative method
  F_new = feval(f,x(ii)+dx,y_approx);
  y_approx = y(:,ii) + dx/2*(F_now + F_new);
  % I stop after one iteration and accept y_approx as y(x+dx).
  y(:,ii+1) = y_approx;
  x(ii+1) = x_0 + ii*dx;
end