% This code applies the most accurate two-step explicit
% multi-step 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 = 1;
Nvec = [10];

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 2-step method
    v(n+1) = -4*v(n) + 5*v(n-1) + dt*(4*f + 2*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);
plot(t,1./(1+t));
xlabel('t');
axis([0,Tmax,-20,20]);