% This code applies the third-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 = 200;
Nvec = [256];


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;
  
  
  % Store old f value (currently in f)
  f0 = f;
  
  % Calculate new f
  f = -v(2)^2;
  
  % Update using Adams-Bashforth 2 (AB2)
  v(3) = v(2) + dt*(1.5*f - 0.5*f0);
  
  % Start iterative loop
  for n = 3:N,

    % Store old f valuew
    f1 = f0;
    f0 = f;
    
    % Calculate new f
    f = -v(n)^2;

    % Update using Adams-Bashforth 2 (AB2)
    v(n+1) = v(n) + dt*(23/12*f - 16/12*f0 + 5/12*f1);

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

end

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