%
% Main program ode3x3 computes the numerical
% solution to the 3 x 3 ODE system by eight 
% integrators
%
% Step through ten integrators
  for int=1:10
%
% Integration parameters
  [n,nout,nsteps,t0,tf,abserr,relerr]=intpar;
%
% Initial condition vector
  [u0]=inital(n,t0);
%
% Output interval  
  dt=tf-t0;
%
% Compute solution at nout output points  
  for j=1:nout
%   
%   Print current solution   
    [out]=fprint(int,n,t0,u0);
%
%   Fixed step modified Euler integrator
    if int == 1
      [u]=euler2a(n,t0,tf,u0,nsteps);
    end
%
%   Variable step modified Euler integrator
    if int == 2
      [u]=euler2b(n,t0,tf,u0,nsteps,abserr,relerr);
    end
%
%   Fixed step classical fourth order RK integrator
    if int == 3
      [u]=rkc4a(n,t0,tf,u0,nsteps);
    end
%
%   Variable step classical fourth order RK integrator
    if int == 4
      [u]=rkc4b(n,t0,tf,u0,nsteps,abserr,relerr);
    end
%
%   Fixed step RK Fehlberg (RKF45) integrator
    if int == 5
      [u]=rkf45a(n,t0,tf,u0,nsteps);
    end
%
%   Variable step RK Fehlberg (RKF45) integrator
    if int == 6
      [u]=rkf45b(n,t0,tf,u0,nsteps,abserr,relerr);
    end
%
%   Fixed step RK (2,3) (embedded pair) integrator
    if int == 7
      [u]=ode23a(n,t0,tf,u0,nsteps);
    end
%
%   Variable step RK (2,3) (embedded pair) integrator
    if int == 8
      [u]=ode23b(n,t0,tf,u0,nsteps,abserr,relerr);
    end
%
%   Fixed step RK (4,5) (embedded pair) integrator
    if int == 9
      [u]=ode45a(n,t0,tf,u0,nsteps);
    end
%
%   Variable step RK (4,5) (embedded pair) integrator
    if int == 10
      [u]=ode45b(n,t0,tf,u0,nsteps,abserr,relerr);
    end
%
%   Advance solution
    t0=tf;
    tf=tf+dt;
    u0=u;
%
% Next output
  end
%
% Next integrator
  end
%
% End of ode3x3
