% This code applies the fourth-order Runge-Kutta 
% 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 = [250,500,1000];

e = zeros(3,100+1);

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

  v = zeros(1,N+1);
  t = linspace(0,Tmax,N+1);
  v(1) = 1.0;
  
  % Start iterative loop
  for n = 1:N,

    % 1st-stage
    a = -dt*v(n)^2;
    va = v(n) + 0.5*a;
    
    % 2nd-stage
    b = -dt*va^2;
    vb = v(n) + 0.5*b;
    
    % 3rd-stage
    c = -dt*vb^2;
    vc = v(n) + c;
    
    % 4th-stage
    d = -dt*vc^2;
    v(n+1) = v(n) + (1/6)*(a + 2*b + 2*c + d);

  end
  
  e(nnn,1:N+1) = abs(v - 1./(1+t));
  
end

% Normalize e by finest e
for nnn = 1:3,
  enorm = e(nnn,1:Nvec(nnn)+1)./e(3,1:(2^(3-nnn)):Nvec(3)+1);
  plot(linspace(0,Tmax,Nvec(nnn)+1),enorm);hold on;
end
xlabel('time');
ylabel('Error/Fine Error');
grid on;