clear all;
close all;

% Set time length of integration, and number of steps
Tmax = 10;
Nvec = [250,500,1000];

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,

    % Calculate rhs
    f = -v(n)^2;

    % Update using Forward Euler
    v(n+1) = v(n-1) + 2*dt*f;

  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;