function ode_check_error
%ODE_TEST  Run non-adaptive ODE solvers of different orders.
%   ODE_TEST compares the orders of accuracy of several explicit Runge-Kutta
%   methods. The non-adaptive ODE solvers are tested on a problem used in
%   E. Hairer, S.P. Norsett, and G. Wanner, Solving Ordinary Differential
%   Equations I, Nonstiff Problems, Springer-Verlag, 1987. 
%   The errors of numerical solutions obtained with several step sizes is
%   plotted against the step size used. For a given solver, the slope of that 
%   line corresponds to the order of the integration method used.  
%

solvers = { 'ode1', 'ode2', 'ode3', 'ode4', 'ode5' };
Nsteps  = [ 10 20 40 80 160 320 640 1280 ];

x0   = 0 ;
y0   = [1;1;1;1] ;
xend = 1 ;
h    = (xend-x0)./Nsteps ;

err = zeros(length(solvers),length(h));

for i = 1:length(solvers)
  solver = solvers{i};
  for j = 1:length(Nsteps)
    N = Nsteps(j);
    x = linspace(x0,xend,N+1);
    y = feval(solver,@f,x,y0);
    err(i,j) = max(max(abs(y - yexact(x))));
  end
end  

for i = 1:length(solvers)
  disp( ' ' )
  disp( ['estimate order for scheme = ' solvers{i}] );
  disp( sprintf( 'Error = %-15g', err(i,1) ) ) ;
  for j = 2:length(Nsteps)
    order = log( err(i,j-1)/err(i,j) ) / log( h(j-1)/h(j) ) ;
    disp( sprintf( 'Error = %-15g Order = %-15g', err(i,j), order) ) ;
  end
end  

% -------------------------
function yp = f(x,y)
yp = [  2*x*y(4)*y(1)
       10*x*y(4)*y(1)^5
        2*x*y(4)
       -2*x*(y(3)-1) ];

% -------------------------
function y = yexact(x)
x2 = x.^2;
y  = zeros(length(x),4);
y(:,1) = exp(sin(x2));
y(:,2) = exp(5*sin(x2));
y(:,3) = sin(x2)+1;
y(:,4) = cos(x2);

