% This code applies the second-order Adams-Bashforth
% method to solve du/dt = -u^2 with u(0) = 1.  Note: the
% exact solution is u = 1/(1+t).

clear all;
close all;

% Set time length of integration, and number of steps
Tmax = 10;
Nvec = [25,50,100];

for nnn = 1:length(Nvec),
  N = Nvec(nnn);
  dt = Tmax/N;

  v(1) = 1.0;
  
  % Use forward Euler for first step
  f = -v(1)^2;
  v(2) = v(1) + dt*f;
  
  % Start iterative loop
  for n = 2:N,

    % Store old f value (currently in f)
    f0 = f;
    
    % Calculate new f
    f = -v(n)^2;

    % Update using Adams-Bashforth 2 (AB2)
    v(n+1) = v(n) + dt*(1.5*f - 0.5*f0);

  end
  
  % Plot results
  t = linspace(0,Tmax,N+1);
  subplot(211);
  plot(t,v,'*');hold on;
  ylabel('v');

  subplot(212);
  e = abs(v - 1./(1+t));
  plot(t,e,'*');hold on;
  ylabel('e');
  
end

% Plot exact solution
t = linspace(0,Tmax,1001);
subplot(211);
plot(t,1./(1+t));

subplot(212);
grid on;