From octave-sources-request at bevo dot che dot wisc dot edu Sun Aug  4 22:00:11 2002
Subject: ode solver source updates
From: Marc Compere <CompereM at asme dot org>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Sun, 4 Aug 2002 17:32:21 -0500


--=-+D8QRHJhVFcPnLLnycDs
Content-Type: text/plain
Content-Transfer-Encoding: 7bit

There are 7 files included in text-format below.  Each file is separated
by dashes.
---------------------------------------------------------
Readme.txt:
----------

This is the Readme.txt file for Octave m-file ordinary differential
equation (ODE) solvers, version 1.15.  They also work in Matlab.

This directory contains 8 files that provide an Octave user with
several options for numerically integrating ODE.
There are 3 fixed-step Runge-Kutta algorithms and
3 variable step Runge-Kutta-Fehlberg algorithms along with
a Dormand-Prince 4(5) pair used by default in ode45.m (as of
version 1.09.)
All are explicit RK formulas that work well with nonstiff or mildly
stiff problems.
All contain their own documentation accessable at the Octave (or
Matlab) command prompt by typing 'help ode45' or 'help rk8fixed'
or whichever solver you need help with.


----------------------------------------------------------------------
The archive ode_v1.15.tar.gz contains 6 explicit single-step
Runge-Kutta ODE solvers along with 2 files demonstrating example
uses of each solver:

   - ode23.m    : variable step, 2nd-3rd order
   - ode45.m    : variable step, 4th-5th order
   - ode78.m    : variable step, 7th-8th order

   - rk2fixed.m : fixed step, 2nd order
   - rk4fixed.m : fixed step, 4th order
   - rk8fixed.m : fixed step, 8th order

   - pendulum.m : a sample m-file script that runs all solvers
   - penddot.m  : derivative function file, returning dy/dt for a simple pendulum
   - Readme.txt : this file

----------------------------------------------------------------------
Steps for testing these ode solvers in Octave (version 2.0.15 or better),
from a unix shell:

(1) unzip and untar the archive:
	gunzip ode_v1.15.tar.gz
	tar xvf ode_v1.15.tar
(2) change directories into the newly created directory and execute Octave:
	cd ode_v1.15
	octave
(3) run the sample pemdulum script from within Octave.
	pendulum

This script sequentially executes all 6 m-file integrators and, if plotting
capability is setup properly, plots the output from all integrators.

----------------------------------------------------------------------
I've made an effort to make these portable to most Octave installations
as well as for use in most Matlab versions.  These work with no modification
in Matlab v5.2 & v5.3.  If you want to use these in Matlab, however,
you'll do yourself a favor by renaming ode45.m and ode23.m to something else,
like ode45_octave.m and ode23_octave.m.  This is because Matlab already has
two integrators named ode45 and ode23.
Don't forget to change the function names at the top of ode45_octave.m and
ode23_octave.m as well.

Some effort has been made to create ode45 and ode23 with similar
argument structures to their Matlab counterparts.  Only the most basic function
calls to ode45 and ode23 match up in both Matlab's integrators and these.
Feel free to change them as you see fit.  You are welcome to mail me for useful
change suggestions.

----------------------------------------------------------------------
Basic differences between the integrators:

In general, the higher the integration order, the smaller the local truncation
error is at each time step.  Small local truncation errors result in larger
integration steps.  This is demonstrated by ode78 generating far fewer
steps than ode23 for solving the same problem over the same time interval with
the same error criterion.
The cost of the higher order integrators is the number of function evaluations
required at each step.  This results in longer execution times for each integration
step however, if the higher order integrator takes much larger steps, the
total computational cost is less.
The tradeoff between solving a problem with an integrator that takes fewer steps
versus using one that takes more time for each step will vary with each problem.
Factors such as numerical stiffness (or lack of) and the number of discontinuities
present in the ODE's will somteimes cause ode23 to be more effective than ode78,
& vice-versa.

Integrator costs in right-hand-side (RHS) evaluations:
   - ode23.m    : requires 3 RHS function evaluations per step
   - ode45.m    : requires 6 RHS function evaluations per step (+1 for the Dormand-Prince pair)
   - ode78.m    : requires 13 RHS function evaluations per step

   - rk2fixed.m : requires 2 RHS function evaluations per step
   - rk4fixed.m : requires 4 RHS function evaluations per step
   - rk8fixed.m : requires 13 RHS function evaluations per step

If you are interested in experimenting with pendulum.m, try turning on
the 'trace' variable for screen output.  This increases execution time
but provides a way to monitor the problem during the integration.

----------------------------------------------------------------------
Two of the original files rk4fixed.m and ode78.m came from other
people.

rk4fixed.m was written by:
     Dr. Raul Longoria, Dept. of Mechanical Engineering,
     The Univ. of Texas at Austin and

ode78.m was originally written by:
     Dr. Howard Wilson & 
     Daljeet Singh, Dept. Of Electrical Engineering,
     The University of Alabama

The other files were created from the structure of these 
along with coefficients from standard numerical methods books.
Dr. Longoria has given permission to redistribute rk4fixed.m.

ode78.m was originally found at:
ftp://ftp.mathworks.com/pub/contrib/v4/diffeq/ode78.m.
& I am redistributing a modified version with Dr. Wilson's
permission.  

Numerous applications in ordinary and partial differential equations
can be found in Dr. Wilson's text:

     Howard Wilson and Louis Turcotte, 'Advanced Mathematics and 
     Mechanics Applications Using MATLAB', 2nd Ed, CRC Press, 1997

----------------------------------------------------------------------

Bugs or changes or comments should be directed to the email address
below.  The latest and greatest versions will generally _NOT_ be found
on the octave-source list, but rather will be maintained  at:
	http://marc.me.utexas.edu/tmp/octave_ode_solvers/
and hopefully sometime soon at:
	http://octave.sourceforge.net/


Marc Compere
CompereM at asme dot org
created : 06 October 1999
modified: 04 August 2002


---------------------------------------------------------
ode23.m:
-------

function [tout,xout] = ode23(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% ode23.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% ode23.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% ode23 (v1.14) Integrates a system of ordinary differential equations using
% 2nd & 3rd order Runge-Kutta formulas.  This particular 3rd-order method reduces
% to Simpson's 1/3 rule and uses the 3rd order estimate for xout.
%
% 3rd-order accurate RK methods have local and global errors of O(h^4) and O(h^3),
% respectively and yield exact results when the solution is a cubic.
%
% The order of the RK method is the order of the local *truncation* error, d,
% which is the principle error term in the portion of the Taylor series
% expansion that gets dropped, or intentionally truncated.  This is different
% from the local error which is the difference between the estimated solution
% and the actual, or true solution.  The local error is used in stepsize
% selection and may be approximated by the difference between two estimates of
% different order, l(h) = x_(O(h+1)) - x_(O(h)).  With this definition, the
% local error will be as large as the error in the lower order method.
% The local truncation error is within the group of terms that gets multipled
% by h when solving for a solution from the general RK method.  Therefore, the
% order-p solution created by the RK method will be roughly accurate to O(h^(p+1))
% since the local truncation error shows up in the solution as h*d, which is
% h times an O(h^(p)) term, or rather O(h^(p+1)).
% Summary:   For an order-p accurate RK method,
%            - the local truncation error is O(h^p)
%            - the local error used for stepsize adjustment and that
%              is actually realized in a solution is O(h^(p+1))
%
% This requires 3 function evaluations per integration step.
%
% Relevant discussion on step size choice can be found on pp.90,91 in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential Equations
% and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% The error estimate formula and slopes are from
% Numerical Methods for Engineers, 2nd Ed., Chapra & Canale, McGraw-Hill, 1985
%
% Usage:
%         [tout, xout] = ode23(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)
%
% INPUT:
% FUN   - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where FUN = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - Initial value COLUMN-vector.
% ode_fcn_format - this specifies if the user-defined ode function is in
%         the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%         or:           xprime = fun(x,t)   (ode_fcn_format=1)
%         Matlab's solvers comply with ode_fcn_format=0 while
%         Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% tol   - The desired accuracy. (optional, default: tol = 1.e-6).
% trace - If nonzero, each step is printed. (optional, default: trace = 0).
% count - if nonzero, variable 'rhs_counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'rhs_counter' is incremented in here, not in the state-dot file
%         simply make 'rhs_counter' global in the file that calls ode23
% hmax  - limit the maximum stepsize to be less than or equal to hmax
%
% OUTPUT:
% tout  - Returned integration time points (column-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 27 June 2001

if nargin < 8, hmax = (tspan(2) - tspan(1))/2.5; end
if nargin < 7, count = 0; end
if nargin < 6, trace = 0; end
if nargin < 5, tol = 1.e-3; end
if nargin < 4, ode_fcn_format = 0; end

pow = 1/4; % see p.91 in the Ascher & Petzold reference for more infomation.

% The 2(3) coefficients:
 a(1,1)=0;
 a(2,1)=1/2;
 a(3,1)=-1; a(3,2)=2;
 % 2nd order b-coefficients
 b2(1)=0; b2(2)=1;
 % 5th order b-coefficients
 b3(1)=1/6; b3(2)=2/3; b3(3)=1/6;
 for i=1:3
  c(i)=sum(a(i,:));
 end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmin = (tfinal - t)/1e12;
h = (tfinal - t)/200; % initial guess at a step size
x = x0(:);            % this always creates a column vector, x
tout = t;             % first output time
xout = x.';           % first output solution

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0; end
end % if count

if trace
 clc, t, h, x
end

% The main loop
   while (t < tfinal) & (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % compute the slopes
      if (ode_fcn_format==0), % (default)
         k(:,1)=feval(FUN,t,x);
         k(:,2)=feval(FUN,t+c(2)*h,x+h*(a(2,1)*k(:,1)));
         k(:,3)=feval(FUN,t+c(3)*h,x+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)));
      else, % ode_fcn_format==1
         k(:,1)=feval(FUN,x,t);
         k(:,2)=feval(FUN,x+h*(a(2,1)*k(:,1)),t+c(2)*h);
         k(:,3)=feval(FUN,x+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)),t+c(3)*h);
      end % if (ode_fcn_format==1)

      % increment rhs_counter
      if count==1, rhs_counter = rhs_counter + 3; end

      % compute the 2nd order estimate
      x2=x + h*b2(2)*k(:,2);
      % compute the 3rd order estimate
      x3=x + h*(b3(1)*k(:,1) + b3(2)*k(:,2) + b3(3)*k(:,3));

      % estimate the local truncation error
      gamma1 = x3 - x2;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');
      tau = tol*max(norm(x,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         x = x3;    % <-- using the higher order estimate is called 'local extrapolation'
         tout = [tout; t];
         xout = [xout; x.'];
      end
      if trace
         home, t, h, x
      end

      % Update the step size
      if delta == 0.0
       delta = 1e-16;
      end
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end;

   if (t < tfinal)
      disp('Step size grew too small.')
      t, h, x
   end


---------------------------------------------------------
ode45.m:
-------

function [tout,xout,Nsteps_acc,Nsteps_rej] = ode45(FUN,tspan,x0,pair,ode_fcn_format,tol,trace,count,hmax,N_est_acc_steps)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% ode45.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% ode45.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% ode45 (v1.15) integrates a system of ordinary differential equations using
% 4th & 5th order embedded formulas from Dormand & Prince or Fehlberg.
%
% The Fehlberg 4(5) pair is established and works well, however, the
% Dormand-Prince 4(5) pair minimizes the local truncation error in the
% 5th-order estimate which is what is used to step forward (local extrapolation.)
% Generally it produces more accurate results and costs roughly the same
% computationally.  The Dormand-Prince pair is the default.
%
% This is a 4th-order accurate integrator therefore the local error normally
% expected is O(h^5).  However, because this particular implementation
% uses the 5th-order estimate for xout (i.e. local extrapolation) moving
% forward with the 5th-order estimate should yield local error of O(h^6).
%
% The order of the RK method is the order of the local *truncation* error, d,
% which is the principle error term in the portion of the Taylor series
% expansion that gets dropped, or intentionally truncated.  This is different
% from the local error which is the difference between the estimated solution
% and the actual, or true solution.  The local error is used in stepsize
% selection and may be approximated by the difference between two estimates of
% different order, l(h) = x_(O(h+1)) - x_(O(h)).  With this definition, the
% local error will be as large as the error in the lower order method.
% The local truncation error is within the group of terms that gets multipled
% by h when solving for a solution from the general RK method.  Therefore, the
% order-p solution created by the RK method will be roughly accurate to O(h^(p+1))
% since the local truncation error shows up in the solution as h*d, which is
% h times an O(h^(p)) term, or rather O(h^(p+1)).
% Summary:   For an order-p accurate RK method,
%            - the local truncation error is O(h^p)
%            - the local error used for stepsize adjustment and that
%              is actually realized in a solution is O(h^(p+1))
%
% This requires 6 function evaluations per integration step.
%
% Both the Dormand-Prince and Fehlberg 4(5) coefficients are from a tableu in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential Equations
% and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% The error estimate formula and slopes are from
% Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle, McGraw-Hill, 1985
%
% Usage:
%         [tout, xout] = ode45(FUN,tspan,x0,pair,ode_fcn_format,tol,trace,count,hmax)
%
% INPUTS:
%    FUN   - String containing name of user-supplied problem description.
%            Call: xprime = fun(t,x) where FUN = 'fun'.
%            t      - Time (scalar).
%            x      - Solution column-vector.
%            xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
%    tspan - [ tstart, tfinal ]
%    x0    - Initial value COLUMN-vector.
%    pair  - flag specifying which integrator coefficients to use:
%               0 --> use Dormand-Prince 4(5) pair (default)
%               1 --> use Fehlberg pair 4(5) pair
%    ode_fcn_format - this specifies if the user-defined ode function is in
%            the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%            or:           xprime = fun(x,t)   (ode_fcn_format=1)
%            Matlab's solvers comply with ode_fcn_format=0 while
%            Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
%    tol   - The desired accuracy. (optional, default: tol = 1.e-6).
%    trace - If nonzero, each step is printed. (optional, default: trace = 0).
%    count - if nonzero, variable 'rhs_counter' is initalized, made global
%            and counts the number of state-dot function evaluations
%            'rhs_counter' is incremented in here, not in the state-dot file
%            simply make 'rhs_counter' global in the file that calls ode45
%    hmax  - limit the maximum stepsize to be less than or equal to hmax
%    N_est_acc_steps - an estimate of how many accepted steps there may be;
%                      used to preallocate memory for the [tout,xout] solutions
%
% OUTPUTS:
%    tout  - array of discretized times points (an Nsteps_acc by 1 column-vector).
%    xout  - solution, one solution column-vector per tout-value (an Nsteps_acc by Nstates matrix)
%    Nsteps_acc - total number of accepted integration steps + 1
%    Nsteps_rej - total number of rejected integration steps
%
% The result can be displayed by: plot(tout, xout).
%
% The following relationship should hold after a completed run:
%    rhs_counter == (Nsteps_acc-1+Nsteps_rej)*6+1
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 03 July 2001

if nargin < 10, N_est_acc_steps = (tspan(2)-tspan(1))*1e3; end
if nargin <  9, hmax = (tspan(2) - tspan(1))/2.5; end
if nargin <  8, count = 0; end
if nargin <  7, trace = 0; end
if nargin <  6, tol = 1.e-6; end
if nargin <  5, ode_fcn_format = 0; end
if nargin <  4, pair = 0; end

pow = 1/6; % see p.91 in the Ascher & Petzold reference for more infomation.

if (pair==0),
   % Use the Dormand-Prince 4(5) coefficients:
    a_(1,1)=0;
    a_(2,1)=1/5;
    a_(3,1)=3/40; a_(3,2)=9/40;
    a_(4,1)=44/45; a_(4,2)=-56/15; a_(4,3)=32/9;
    a_(5,1)=19372/6561; a_(5,2)=-25360/2187; a_(5,3)=64448/6561; a_(5,4)=-212/729;
    a_(6,1)=9017/3168; a_(6,2)=-355/33; a_(6,3)=46732/5247; a_(6,4)=49/176; a_(6,5)=-5103/18656;
    a_(7,1)=35/384; a_(7,2)=0; a_(7,3)=500/1113; a_(7,4)=125/192; a_(7,5)=-2187/6784; a_(7,6)=11/84;
    % 4th order b-coefficients
    b4_(1,1)=5179/57600; b4_(2,1)=0; b4_(3,1)=7571/16695; b4_(4,1)=393/640; b4_(5,1)=-92097/339200; b4_(6,1)=187/2100; b4_(7,1)=1/40;
    % 5th order b-coefficients
    b5_(1,1)=35/384; b5_(2,1)=0; b5_(3,1)=500/1113; b5_(4,1)=125/192; b5_(5,1)=-2187/6784; b5_(6,1)=11/84; b5_(7,1)=0;
    for i=1:7,
     c_(i)=sum(a_(i,:));
    end
else, % pair==1 so use the Fehlberg 4(5) coefficients:
    a_(1,1)=0;
    a_(2,1)=1/4;
    a_(3,1)=3/32; a_(3,2)=9/32;
    a_(4,1)=1932/2197; a_(4,2)=-7200/2197; a_(4,3)=7296/2197;
    a_(5,1)=439/216; a_(5,2)=-8; a_(5,3)=3680/513; a_(5,4)=-845/4104;
    a_(6,1)=-8/27; a_(6,2)=2; a_(6,3)=-3544/2565; a_(6,4)=1859/4104; a_(6,5)=-11/40;
    % 4th order b-coefficients (guaranteed to be a column vector)
    b4_(1,1)=25/216; b4_(2,1)=0; b4_(3,1)=1408/2565; b4_(4,1)=2197/4104; b4_(5,1)=-1/5;
    % 5th order b-coefficients (also guaranteed to be a column vector)
    b5_(1,1)=16/135; b5_(2,1)=0; b5_(3,1)=6656/12825; b5_(4,1)=28561/56430; b5_(5,1)=-9/50; b5_(6,1)=2/55;
    for i=1:6,
     c_(i)=sum(a_(i,:));
    end
end % if (pair==0)


% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmin = (tfinal - t)/1e20;
h = (tfinal - t)/100; % initial step size guess
x = x0(:);            % ensure x is a column vector

Nstates = size(x,1);
tout = zeros(N_est_acc_steps,1);       % preallocating memory is not only faster but, for long
xout = zeros(N_est_acc_steps,Nstates); % simulations, prevents disappointing surprises later

Nsteps_rej = 0;
Nsteps_acc = 1;
tout(Nsteps_acc) = t;     % first output time
xout(Nsteps_acc,:) = x.'; % first output solution = IC's

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0; end
end % if count

if trace
 clc, t, h %, x
end

if (pair==0),

   k_ = x*zeros(1,7);  % k_ needs to be initialized as an Nx7 matrix where N=number of rows in x
                       % (just for speed so octave doesn't need to allocate more memory at each stage value)

   % Compute the first stage prior to the main loop.  This is part of the Dormand-Prince pair caveat.
   % Normally, during the main loop the last stage for x(k) is the first stage for computing x(k+1).
   % So, the very first integration step requires 7 function evaluations, then each subsequent step
   % 6 function evaluations because the first stage is simply assigned from the last step's last stage.
   % note: you can see this by the last element in c_ is 1.0, thus t+c_(7)*h = t+h, ergo, the next step.
   if (ode_fcn_format==0), % (default)
      k_(:,1)=feval(FUN,t,x); % first stage
   else, % ode_fcn_format==1
      k_(:,1)=feval(FUN,x,t);
   end % if (ode_fcn_format==1)

   % increment rhs_counter
   if (count==1), rhs_counter = rhs_counter + 1; end

   % The main loop using Dormand-Prince 4(5) pair
   while (t < tfinal) & (h >= hmin),
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes by computing the k(:,j+1)'th column based on the previous k(:,1:j) columns
      % notes: k_ needs to end up as an Nxs, a_ is 7x6, which is s by (s-1),
      %        s is the number of intermediate RK stages on [t (t+h)] (Dormand-Prince has s=7 stages)
      if (ode_fcn_format==0), % (default)
         for j = 1:6,
            k_(:,j+1) = feval(FUN, t+c_(j+1)*h, x+h*k_(:,1:j)*a_(j+1,1:j)');
         end
      else, % ode_fcn_format==1
         for j = 1:6,
            k_(:,j+1) = feval(FUN, x+h*k_(:,1:j)*a_(j+1,1:j)', t+c_(j+1)*h);
         end
      end % if (ode_fcn_format==1)

      % increment rhs_counter
      if (count==1), rhs_counter = rhs_counter + 6; end

      % compute the 4th order estimate
      x4=x + h* (k_*b4_); % k_ is Nxs (or Nx7) and b4_ is a 7x1

      % compute the 5th order estimate
      x5=x + h*(k_*b5_); % k_ is Nxs (or Nx7) and b5_ is a 7x1

      % estimate the local truncation error
      gamma1 = x5 - x4;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');       % actual error
      tau = tol*max(norm(x,'inf'),1.0); % allowable error

      % Update the solution only if the error is acceptable
      if (delta<=tau),

         t = t + h;
         x = x5;    % <-- using the higher order estimate is called 'local extrapolation'
         Nsteps_acc = Nsteps_acc + 1;
         tout(Nsteps_acc) = t;
         xout(Nsteps_acc,:) = x.';

         % Assign the last stage for x(k) as the first stage for computing x(k+1).
         % This is part of the Dormand-Prince pair caveat.
         % k_(:,7) has already been computed, so use it instead of recomputing it
         % again as k_(:,1) during the next step.
         k_(:,1)=k_(:,7);

      else, % unacceptable integration step

         Nsteps_rej = Nsteps_rej + 1;

      end % if (delta<=tau)

      if trace
         home, t, h, Nsteps_acc, Nsteps_rej
      end % if trace

      % Update the step size
      if (delta==0.0),
       delta = 1e-16;
      end % if (delta==0.0)
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end % while (t < tfinal) & (h >= hmin)

else, % pair==1 --> use RK-Fehlberg 4(5) pair

   k_ = x*zeros(1,6);  % k_ needs to be initialized as an Nx6 matrix where N=number of rows in x
                       % (just for speed so octave doesn't need to allocate more memory at each stage value)
 
   % The main loop using Dormand-Prince 4(5) pair
   while (t < tfinal) & (h >= hmin),
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes by computing the k(:,j+1)'th column based on the previous k(:,1:j) columns
      % notes: k_ needs to end up as an Nx6, a_ is 6x5, which is s by (s-1),  (RK-Fehlberg has s=6 stages)
      %        s is the number of intermediate RK stages on [t (t+h)]
      if (ode_fcn_format==0), % (default)
         k_(:,1)=feval(FUN,t,x); % first stage
      else, % ode_fcn_format==1
         k_(:,1)=feval(FUN,x,t);
      end % if (ode_fcn_format==1)

      if (ode_fcn_format==0), % (default)
         for j = 1:5,
            k_(:,j+1) = feval(FUN, t+c_(j+1)*h, x+h*k_(:,1:j)*a_(j+1,1:j)');
         end
      else, % ode_fcn_format==1
         for j = 1:5,
            k_(:,j+1) = feval(FUN, x+h*k_(:,1:j)*a_(j+1,1:j)', t+c_(j+1)*h);
         end
      end % if (ode_fcn_format==1)

      % increment rhs_counter
      if (count==1), rhs_counter = rhs_counter + 6; end

      % compute the 4th order estimate
      x4=x + h* (k_(:,1:5)*b4_); % k_(:,1:5) is an Nx5 and b4_ is a 5x1

      % compute the 5th order estimate
      x5=x + h*(k_*b5_); % k_ is the same as k_(:,1:6) and is an Nx6 and b5_ is a 6x1

      % estimate the local truncation error
      gamma1 = x5 - x4;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');       % actual error
      tau = tol*max(norm(x,'inf'),1.0); % allowable error

      % Update the solution only if the error is acceptable
      if (delta<=tau),

         t = t + h;
         x = x5;    % <-- using the higher order estimate is called 'local extrapolation'
         Nsteps_acc = Nsteps_acc + 1;
         tout(Nsteps_acc) = t;
         xout(Nsteps_acc,:) = x.';

      else, % unacceptable integration step

         Nsteps_rej = Nsteps_rej + 1;

      end % if (delta<=tau)

      if trace
         home, t, h, Nsteps_acc, Nsteps_rej
      end % if trace

      % Update the step size
      if (delta==0.0),
       delta = 1e-16;
      end % if (delta==0.0)
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end % while (t < tfinal) & (h >= hmin)

end % if (pair==0),


% if necessary, trim outputs
if (Nsteps_acc<N_est_acc_steps),
   tout = tout(1:Nsteps_acc);
   xout = xout(1:Nsteps_acc,:);
end % if (Nsteps_acc<N_est_acc_steps)


if (t < tfinal),
   disp('Step size grew too small.')
   t, h, Nsteps_acc, Nsteps_rej
end


---------------------------------------------------------
ode78.m:
-------

function [tout,xout] = ode78(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% ode78.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% ode78.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% ode78 (v1.14) Integrates a system of ordinary differential equations using
% 7th order formulas.
%
% This is a 7th-order accurate integrator therefore the local error normally
% expected is O(h^8).  However, because this particular implementation
% uses the 8th-order estimate for xout (i.e. local extrapolation) moving
% forward with the 8th-order estimate will yield errors on the order of O(h^9).
%
% The order of the RK method is the order of the local *truncation* error, d,
% which is the principle error term in the portion of the Taylor series
% expansion that gets dropped, or intentionally truncated.  This is different
% from the local error which is the difference between the estimated solution
% and the actual, or true solution.  The local error is used in stepsize
% selection and may be approximated by the difference between two estimates of
% different order, l(h) = x_(O(h+1)) - x_(O(h)).  With this definition, the
% local error will be as large as the error in the lower order method.
% The local truncation error is within the group of terms that gets multipled
% by h when solving for a solution from the general RK method.  Therefore, the
% order-p solution created by the RK method will be roughly accurate to O(h^(p+1))
% since the local truncation error shows up in the solution as h*d, which is
% h times an O(h^(p)) term, or rather O(h^(p+1)).
% Summary:   For an order-p accurate RK method,
%            - the local truncation error is O(h^p)
%            - the local error used for stepsize adjustment and that
%              is actually realized in a solution is O(h^(p+1))
%
% This requires 13 function evaluations per integration step.
%
% Relevant discussion on step size choice can be found on pp.90,91 in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential Equations
% and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% More may be found in the original author's text containing numerous
% applications on ordinary and partial differential equations using Matlab:
%
%     Howard Wilson and Louis Turcotte, 'Advanced Mathematics and 
%     Mechanics Applications Using MATLAB', 2nd Ed, CRC Press, 1997
%
%
% [tout, xout] = ode78(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)
%
% INPUT:
% FUN   - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where FUN = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - Initial value COLUMN-vector.
% ode_fcn_format - this specifies if the user-defined ode function is in
%         the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%         or:           xprime = fun(x,t)   (ode_fcn_format=1)
%         Matlab's solvers comply with ode_fcn_format=0 while
%         Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% tol   - The desired accuracy. (optional, default: tol = 1.e-6).
% trace - If nonzero, each step is printed. (optional, default: trace = 0).
% count - if nonzero, variable 'rhs_counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'rhs_counter' is incremented in here, not in the state-dot file
%         simply make 'rhs_counter' global in the file that calls ode78
% hmax  - limit the maximum stepsize to be less than or equal to hmax
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).

%   Daljeet Singh & Howard Wilson
%   Dept. Of Electrical Engg., The University of Alabama.
%   11-24-1988.
%
% modified by:
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 19 May 2001


% The Fehlberg coefficients:
alpha_ = [ 2./27., 1/9, 1/6, 5/12, 0.5, 5/6, 1/6, 2/3, 1/3, 1, 0, 1 ]';
beta_  = [ 2/27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/36, 1/12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/24, 0, 1/8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          5/12, 0, -25/16, 25/16, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          0.05, 0, 0, 0.25, 0.2, 0, 0, 0, 0, 0, 0, 0, 0 ;
          -25/108, 0, 0, 125/108, -65/27, 125/54, 0, 0, 0, 0, 0, 0, 0 ;
          31/300, 0, 0, 0, 61/225, -2/9, 13/900, 0, 0, 0, 0, 0, 0 ;
          2, 0, 0, -53/6, 704/45, -107/9, 67/90, 3, 0, 0, 0, 0, 0 ;
          -91/108, 0, 0, 23/108, -976/135, 311/54, -19/60, 17/6, -1/12, 0, 0, 0, 0 ;
          2383/4100, 0, 0, -341/164, 4496/1025, -301/82, 2133/4100, 45/82, 45/164, 18/41, 0, 0, 0 ;
          3/205, 0, 0, 0, 0, -6/41, -3/205, -3/41, 3/41, 6/41, 0, 0, 0 ;
          -1777/4100, 0, 0, -341/164, 4496/1025, -289/82, 2193/4100, 51/82, 33/164, 12/41, 0, 1, 0 ]';
chi_  = [ 0, 0, 0, 0, 0, 34/105, 9/35, 9/35, 9/280, 9/280, 0, 41/840, 41/840]';
psi_  = [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1 ]';

pow = 1/8; % see p.91 in the Ascher & Petzold reference for more infomation.


if nargin < 8, hmax = (tspan(2) - tspan(1))/2.5; end
if nargin < 7, count = 0; end
if nargin < 6, trace = 0; end
if nargin < 5, tol = 1.e-6; end
if nargin < 4, ode_fcn_format = 0; end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
% the following step parameters are used in ODE45
% hmax = (tfinal - t)/5;
% hmin = (tfinal - t)/20000;
% h = (tfinal - t)/100;
% The following parameters were taken because the integrator has
% higher order than ODE45. This choice is somewhat subjective.
%hmin = (tfinal - t)/10000;
hmin = (tfinal - t)/1e20;
h = (tfinal - t)/50;
x = x0(:);          % the '(:)' ensures x is initialized as a column vector
f = x*zeros(1,13);  % f needs to be an Nx13 matrix where N=number of rows in x
tout = t;
xout = x.';
tau = tol * max(norm(x,'inf'), 1);

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

if trace
%  clc, t, h, x
%   clc, t, x
   clc, t
end

% The main loop
   while (t < tfinal) & (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes
      if (ode_fcn_format==0), % (default)
         f(:,1) = feval(FUN,t,x);
         for j = 1: 12,
            f(:,j+1) = feval(FUN, t+alpha_(j)*h, x+h*f*beta_(:,j));
         end
      else, % ode_fcn_format==1
         f(:,1) = feval(FUN,x,t);
         for j = 1: 12,
            f(:,j+1) = feval(FUN, x+h*f*beta_(:,j), t+alpha_(j)*h);
         end
      end % if (ode_fcn_format==1)


      % increment rhs_counter
      if count==1, rhs_counter = rhs_counter + 13; end

      % Truncation error term
      gamma1 = h*41/840*f*psi_;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');
      tau = tol*max(norm(x,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         x = x + h*f*chi_;  % this integrator uses local extrapolation
         tout = [tout; t];
         xout = [xout; x.'];
      end
      if trace
%        home, t, h, x
%        home, t, x
         home, t, h
      end

      % Update the step size
      if delta == 0.0
       delta = 1e-16;
      end
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end;

   if (t < tfinal)
      disp('SINGULARITY LIKELY.')
      t
   end


---------------------------------------------------------
penddot.m:
---------

function xdot=penddot(t,x)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% penddot.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% penddot.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% This is an example derivative function file that works
% with Octave or Matlab.
% The equations represent motion of a simple pendulum with damping.
%
% The plots created by pendulum.m show the angular position and velocity
% trajectories created by each different integrator.
% Position is the trace that reaches steady state at 0(rad)
% because of the gravity term, -m*g*l/2*sin(x(1)).
% Velocity reaches a steady state of zero because of the
% damping term, -b*x(2).
%
% Use ode45 to integrate these ODE's
% like this:
%    [t,x] = ode45('penddot',tspan,IC);
%
% x is the state column vector and meant to be used only within this m-file
% This function is meant to return the derivatives of the state variable
% given the state vector and time.
%
% Structure:	xdot = [ x1dot, x2dot, ..., xNdot ]'
%
% eg.  ml^2*thetadd + b*thetad + m*g*l*sin(theta) = 0
%
%	x(1) = theta
%	x(2) = thetad 	( = x(1)dot )
%
% Convention: the lowest order states are first columnwise
%
% Marc Compere
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 23 October 2001

global m g l b counter index

xdot=[ x(2) , 1/(1/3*m*l^2)*(-b*x(2)-m*g*l/2*sin(x(1)))]';

% remember to return a column vector

end


---------------------------------------------------------
pendulum.m:
----------

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% pendulum.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% pendulum.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% This integrates a set of ordinary differential equations (ODE) using 6 different
% ODE solvers.  The equations represent the dynamics of a simple pendulum.
%
% The integrators ode78.m, ode45.m, ode23.m, rk8fixed.m, rk4fixed.m, and rk2fixed.m
% all produce column vector output similar to Matlab.
%
% All integrators work in octave 2.1.32 and Matlab 5.3 with no modification.
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 23 October 2001

clear

% allow global access to the parameters m, g, l, & b:
global m g l b rhs_counter

m=1;    % (kg)
g=9.81; % (m/s^2)
l=1.;   % (m)
b=0.7;   % ((N-s)/m))

% integrator setup:
trace = 0;   % this is a (1/0) flag that puts output to the screen or not
count = 1;   % this is a (1/0) flag that causes rk4fixed to increment 'rhs_counter' or not
rhs_counter = 0; % 'rhs_counter' is the number of right-hand-side function evaluations, z_dot=f(z)
ode_fcn_format = 0; % 0 chooses Matlab-format right-hand-side function definitions
                    % xdot=f(t,x), or Octave's lsode format, xdot=f(x,t)
t0=0;
tfinal = 5; % (seconds)
tspan = [t0,tfinal];
hmax = 0.1;

% Initial Conditions: theta(t=0)=30(deg) & initially at rest
IC = [ 30*pi/180,0]'; % (rad), (rad/s)
sps = 100;            % sps -> step per second
Nsteps=(tfinal-t0)*sps;    % this creates sps number of integration steps per second
tolerance = 1e-3;

% polite housekeeping
this_user_uses=page_screen_output;
page_screen_output=0;

% for development:
tol=tolerance; pair=0; x0=IC; FUN='penddot';

% Solve the ODE specified in penddot.m using each of the 6 m-file integrators.

t_begin_calcs=cputime;

   disp('Integrating using rk2fixed...')
   [t1,zrk2fixed] = rk2fixed('penddot',tspan,IC,Nsteps,ode_fcn_format,trace,count);     % fixed step integration
   rk2_counter = rhs_counter;
   rhs_counter=0;

t(1)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using rk4fixed...')
   [t2,zrk4fixed] = rk4fixed('penddot',tspan,IC,Nsteps,ode_fcn_format,trace,count);     % fixed step integration
   rk4_counter = rhs_counter;
   rhs_counter=0;

t(2)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using rk8fixed...')
   [t3,zrk8fixed] = rk8fixed('penddot',tspan,IC,Nsteps,ode_fcn_format,trace,count);     % fixed step integration
   rk8_counter = rhs_counter;
   rhs_counter=0;

t(3)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode23...')
   [t4,zode23] = ode23('penddot',tspan,IC,ode_fcn_format,tolerance,trace,count,hmax); % rk23 variable step integration
   ode23_counter = rhs_counter;
   rhs_counter=0;

t(4)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode45 with the Dormand-Prince 4(5) pair...')
   pair=0;
   [t5,zode45dp,Nsteps_acc_ode45dp,Nsteps_rej_ode45dp] = ode45('penddot',tspan,IC,pair,ode_fcn_format,tolerance,trace,count,hmax); % rk45 variable step integration
   ode45dp_counter = rhs_counter;
   rhs_counter=0;

t(5)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode45 with the Fehlberg 4(5) pair...')
   pair=1;
   [t6,zode45rkf,Nsteps_acc_ode45rkf,Nsteps_rej_ode45rkf] = ode45('penddot',tspan,IC,pair,ode_fcn_format,tolerance,trace,count,hmax); % rk45 variable step integration
   ode45rkf_counter = rhs_counter;
   rhs_counter=0;

t(6)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode78...')
   [t7,zode78] = ode78('penddot',tspan,IC,ode_fcn_format,tolerance,trace,count,hmax); % rk78 variable step integration
   ode78_counter = rhs_counter;
   rhs_counter=0;

t(7)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   %disp('Integrating using sdirk...')
   % Note: If you want to use sdirk you have to compile sdirk.oct, then change the function definition in
   %       penddot.m to zdot=penddot(z,t) and uncomment the lines below.
   %       Then if you want the other integrators to still work, set ode_fcn_format=1 above.
   %skip_step=10;
   %h_initial = 1e-3;
   %[t8,zsdirk] = sdirk('penddot',IC,tspan,tolerance,0.1*tolerance,h_initial,trace,skip_step); % 4th order stiff integration
   %sdirk_counter = rhs_counter;
   %rhs_counter=0;

t(8)=cputime-t_begin_calcs;

rhs_fcn_evaluation_summary=[
strcat('rk2_counter      = ',num2str(rk2_counter));
strcat('rk4_counter      = ',num2str(rk4_counter));
strcat('rk8_counter      = ',num2str(rk8_counter));
strcat('ode23_counter    = ',num2str(ode23_counter));
strcat('ode45dp_counter  = ',num2str(ode45dp_counter));
strcat('ode45rkf_counter = ',num2str(ode45rkf_counter));
strcat('ode78_counter    = ',num2str(ode78_counter)) ]
%strcat('sdirk_counter    = ',num2str(sdirk_counter)) ]


disp('Elapsed times for each solver to integrate the state equations for a simple pendulum:')
t


% notes: The Dormand-Prince pair in ode45 produces a time output vector
%        of 18x1.  (size(t5)=18x1)
%        The number of function evaluations is 103 which you can compute
%        from (18-1)*6+1=103.  (18-1) because there were really only 17 new
%        steps computed.  The first step is just the initial conditions.
%        Multiply (18-1) by 6 because each step during the normal main loop
%        requires 6 function evaluations to compute.  Add 1 because the very
%        first step requires 1 function evaluation to start the main loop.


% plot those puppies
figure(1)
clg
if ishold~=1, hold, end
title('Pendulum Position & Velocity')
ylabel('Theta & Theta_dot (rad) & (rad/s)')
xlabel('Time (s)')
plot(t1,zrk2fixed)
plot(t2,zrk4fixed)
plot(t3,zrk8fixed)
plot(t4,zode23)
plot(t5,zode45dp)
plot(t6,zode45rkf)
plot(t7,zode78)
%plot(t8,zsdirk)
if ishold==1, hold, end


% These plots show the angular position and velocity
% trajectories created by each different integrator.
% Position is the trace that reaches steady state at -pi/2.
% Velocity reaches a steady state of zero.


% Return setting(s) to what they were before
page_screen_output=this_user_uses;



---------------------------------------------------------
rk2fixed.m:
----------

function [tout,xout] = rk2fixed(FUN,tspan,x0,Nsteps,ode_fcn_format,trace,count)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% rk2fixed.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% rk2fixed.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% rk2fixed (v1.14) integrates a system of ordinary differential equations using a
% 2nd order Runge-Kutta formula called Ralston's method.
% This choice of 2nd order coefficients provides a minimum bound on truncation error.
% For more, see Ralston & Rabinowitz (1978) or 
% Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle, McGraw-Hill, 1985
%
% 2nd-order accurate RK methods have a local error estimate of O(h^3).
%
% rk2fixed() requires 2 function evaluations per integration step.
%
% Usage:
%         [tout, xout] = rk2fixed(FUN, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% FUN    - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where FUN = 'fun'.
%          t      - Time (scalar).
%          x      - Solution column-vector.
%          xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan  - [ tstart, tfinal ]
% x0     - Initial value COLUMN-vector.
% Nsteps - number of steps used to span [ tstart, tfinal ]
% ode_fcn_format - this specifies if the user-defined ode function is in
%          the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%          or:           xprime = fun(x,t)   (ode_fcn_format=1)
%          Matlab's solvers comply with ode_fcn_format=0 while
%          Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% trace  - If nonzero, each step is printed. (optional, default: trace = 0).
% count  - if nonzero, variable 'rhs_counter' is initalized, made global
%          and counts the number of state-dot function evaluations
%          'rhs_counter' is incremented in here, not in the state-dot file
%          simply make 'rhs_counter' global in the file that calls rk2fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 19 May 2001

if nargin < 7, count = 0; end
if nargin < 6, trace = 0; end
if nargin < 5, Nsteps = 400/(tspan(2)-tspan(1)); end % <-- 400 is a guess for a default,
                                                %  try verifying the solution with rk4fixed
if nargin < 4, ode_fcn_format = 0; end

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/Nsteps;
xout(1,:) = x0';           
tout(1)   = t;
x = x0(:);

if trace
 clc, t, h, x
end

% The main loop
h = (tspan(2)-tspan(1))/Nsteps;

for i=1:Nsteps,
     if (ode_fcn_format==0),
      k1 = feval(FUN,t,x);
      k2 = feval(FUN,t+3/4*h,x+3/4*h*k1);
     else,
      k1 = feval(FUN,x,t);
      k2 = feval(FUN,x+3/4*h*k1,t+3/4*h);
     end % if (ode_fcn_format==0)

     % increment rhs_counter
     if count==1,
      rhs_counter = rhs_counter + 2;
     end % if

     t = t + h;
     x = (x+h*(1/3*k1+2/3*k2));
     tout = [tout; t];
     xout = [xout; x.'];

     if trace,
      home, t, h, x
     end

end


---------------------------------------------------------
rk4fixed.m:
----------

function [tout,xout] = rk4fixed(FUN,tspan,x0,Nsteps,ode_fcn_format,trace,count)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% rk4fixed.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% rk4fixed.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% rk4fixed (v1.14) is a 4th order Runge-Kutta numerical integration routine.
% It requires 4 function evaluations per step.
%
% 4th-order accurate RK methods have a local error estimate of O(h^5).
%
%
% Usage:
%         [tout, xout] = rk4fixed(FUN, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% FUN    - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where FUN = 'fun'.
%          t      - Time or independent variable (scalar).
%          x      - Solution column-vector.
%          xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan  - [ tstart, tfinal ]
% x0     - Initial value COLUMN-vector.
% Nsteps - number of steps used to span [ tstart, tfinal ]
% ode_fcn_format - this specifies if the user-defined ode function is in
%          the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%          or:           xprime = fun(x,t)   (ode_fcn_format=1)
%          Matlab's solvers comply with ode_fcn_format=0 while
%          Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% trace  - If nonzero, each step is printed. (optional, default: trace = 0).
% count  - if nonzero, variable 'rhs_counter' is initalized, made global
%          and counts the number of state-dot function evaluations
%          'rhs_counter' is incremented in here, not in the state-dot file
%          simply make 'rhs_counter' global in the file that calls rk4fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 19 May 2001

if nargin < 7, count = 0; end
if nargin < 6, trace = 0; end
if nargin < 5, Nsteps = 100/(tspan(2)-tspan(1)); end % <-- 100 is a guess for a default,
                                                %  try verifying the solution with rk8fixed
if nargin < 4, ode_fcn_format = 0; end

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/Nsteps;
xout(1,:) = x0';
tout(1) = t;
x = x0(:);
halfh = 0.5*h;

if trace
 clc, t, h, x
end

for i=1:Nsteps,
     if (ode_fcn_format==0),
      RK1 = feval(FUN,t,x);
      thalf = t+halfh;
      xtemp = x+halfh*RK1;
      RK2 = feval(FUN,thalf,xtemp);
      xtemp = x+halfh*RK2;
      RK3 = feval(FUN,thalf,xtemp);
      tfull = t+h;
      xtemp = x+h*RK3;
      RK4 = feval(FUN,tfull,xtemp);
     else,
      RK1 = feval(FUN,x,t);
      thalf = t+halfh;
      xtemp = x+halfh*RK1;
      RK2 = feval(FUN,xtemp,thalf);
      xtemp = x+halfh*RK2;
      RK3 = feval(FUN,xtemp,thalf);
      tfull = t+h;
      xtemp = x+h*RK3;
      RK4 = feval(FUN,xtemp,tfull);
     end % if (ode_fcn_format==0)

     % increment rhs_counter
     if count==1,
      rhs_counter = rhs_counter + 4;
     end % if

     t = t + h;
     x = (x+h/6*(RK1+2.0*(RK2+RK3)+RK4));
     tout = [tout; t];
     xout = [xout; x.'];

     if trace,
      home, t, h, x
     end

end



---------------------------------------------------------
rk8fixed.m:
----------

function [tout,xout] = rk8fixed(FUN,tspan,x0,Nsteps,ode_fcn_format,trace,count)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% rk8fixed.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% rk8fixed.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% rk8fixed (v1.14) is an 8th order Runge-Kutta numerical integration routine.
% It requires 13 function evaluations per step.  This is not the most
% efficient 8th order implementation.  It was just the easiest to put
% together as a variant from ode78.m.
%
% 8th-order accurate RK methods have a local error estimate of O(h^9).
%
%
% Usage:
%         [tout, xout] = rk8fixed(FUN, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% FUN    - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where FUN = 'fun'.
%          t      - Time or independent variable (scalar).
%          x      - Solution column-vector.
%          xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan  - [ tstart, tfinal ]
% x0     - Initial value COLUMN-vector.
% Nsteps - number of steps used to span [ tstart, tfinal ]
% ode_fcn_format - this specifies if the user-defined ode function is in
%          the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%          or:           xprime = fun(x,t)   (ode_fcn_format=1)
%          Matlab's solvers comply with ode_fcn_format=0 while
%          Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% trace  - If nonzero, each step is printed. (optional, default: trace = 0).
% count  - if nonzero, variable 'rhs_counter' is initalized, made global
%          and counts the number of state-dot function evaluations
%          'rhs_counter' is incremented in here, not in the state-dot file
%          simply make 'rhs_counter' global in the file that calls rk4fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 19 May 2001

if nargin < 7, count = 0; end
if nargin < 6, trace = 0; end
if nargin < 5, Nsteps = 50/(tspan(2)-tspan(1)); end % <-- 50 is a guess for a default,
                                                %  try verifying the solution with ode78
if nargin < 4, ode_fcn_format = 0; end

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

alpha_ = [ 2./27., 1/9, 1/6, 5/12, 0.5, 5/6, 1/6, 2/3, 1/3, 1, 0, 1 ]';
beta_  = [ 2/27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/36, 1/12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/24, 0, 1/8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          5/12, 0, -25/16, 25/16, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          0.05, 0, 0, 0.25, 0.2, 0, 0, 0, 0, 0, 0, 0, 0 ;
          -25/108, 0, 0, 125/108, -65/27, 125/54, 0, 0, 0, 0, 0, 0, 0 ;
          31/300, 0, 0, 0, 61/225, -2/9, 13/900, 0, 0, 0, 0, 0, 0 ;
          2, 0, 0, -53/6, 704/45, -107/9, 67/90, 3, 0, 0, 0, 0, 0 ;
          -91/108, 0, 0, 23/108, -976/135, 311/54, -19/60, 17/6, -1/12, 0, 0, 0, 0 ;
          2383/4100, 0, 0, -341/164, 4496/1025, -301/82, 2133/4100, 45/82, 45/164, 18/41, 0, 0, 0 ;
          3/205, 0, 0, 0, 0, -6/41, -3/205, -3/41, 3/41, 6/41, 0, 0, 0 ;
          -1777/4100, 0, 0, -341/164, 4496/1025, -289/82, 2193/4100, 51/82, 33/164, 12/41, 0, 1, 0 ]';
chi_   = [ 0, 0, 0, 0, 0, 34/105, 9/35, 9/35, 9/280, 9/280, 0, 41/840, 41/840]';

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/Nsteps;
xout(1,:) = x0';
tout(1) = t;
x = x0(:);
f = x*zeros(1,13);

if trace
 clc, t, h, x
end

for i=1:Nsteps,

     % Compute the slopes
     if (ode_fcn_format==0),
      f(:,1) = feval(FUN,t,x);
      for j = 1:12
         f(:,j+1) = feval(FUN, t+alpha_(j)*h, x+h*f*beta_(:,j));
      end
     else,
      f(:,1) = feval(FUN,x,t);
      for j = 1:12
         f(:,j+1) = feval(FUN, x+h*f*beta_(:,j), t+alpha_(j)*h);
      end
     end % if (ode_fcn_format==0)

     % increment rhs_counter
     if count==1,
      rhs_counter = rhs_counter + 13;
     end % if

     t = t + h;
     x = x + h*f*chi_;
     tout = [tout; t];
     xout = [xout; x.'];

     if trace,
      home, t, h, x
     end

end


---------------------------------------------------------




-- 
_______________________________________________________
 Marc Compere, Ph.D., The University of Texas at Austin
 CompereM at asme dot org, (512) 826-0729,   <><
 http://marc.me.utexas.edu/
_______________________________________________________

--=-+D8QRHJhVFcPnLLnycDs
Content-Type: text/html; charset=utf-8

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 TRANSITIONAL//EN">
<HTML>
<HEAD>
  <META HTTP-EQUIV="Content-Type" CONTENT="text/html; CHARSET=UTF-8">
  <META NAME="GENERATOR" CONTENT="GtkHTML/1.0.3">
</HEAD>
<BODY>
There are 7 files included in text-format below.&nbsp; Each file is separated by dashes.
<BR>
---------------------------------------------------------
<BR>
<TT>Readme.txt:</TT>
<BR>
----------
<PRE>This is the Readme.txt file for Octave m-file ordinary differential
equation (ODE) solvers, version 1.15.  They also work in Matlab.

This directory contains 8 files that provide an Octave user with
several options for numerically integrating ODE.
There are 3 fixed-step Runge-Kutta algorithms and
3 variable step Runge-Kutta-Fehlberg algorithms along with
a Dormand-Prince 4(5) pair used by default in ode45.m (as of
version 1.09.)
All are explicit RK formulas that work well with nonstiff or mildly
stiff problems.
All contain their own documentation accessable at the Octave (or
Matlab) command prompt by typing 'help ode45' or 'help rk8fixed'
or whichever solver you need help with.


----------------------------------------------------------------------
The archive ode_v1.15.tar.gz contains 6 explicit single-step
Runge-Kutta ODE solvers along with 2 files demonstrating example
uses of each solver:

   - ode23.m    : variable step, 2nd-3rd order
   - ode45.m    : variable step, 4th-5th order
   - ode78.m    : variable step, 7th-8th order

   - rk2fixed.m : fixed step, 2nd order
   - rk4fixed.m : fixed step, 4th order
   - rk8fixed.m : fixed step, 8th order

   - pendulum.m : a sample m-file script that runs all solvers
   - penddot.m  : derivative function file, returning dy/dt for a simple pendulum
   - Readme.txt : this file

----------------------------------------------------------------------
Steps for testing these ode solvers in Octave (version 2.0.15 or better),
from a unix shell:

(1) unzip and untar the archive:
	gunzip ode_v1.15.tar.gz
	tar xvf ode_v1.15.tar
(2) change directories into the newly created directory and execute Octave:
	cd ode_v1.15
	octave
(3) run the sample pemdulum script from within Octave.
	pendulum

This script sequentially executes all 6 m-file integrators and, if plotting
capability is setup properly, plots the output from all integrators.

----------------------------------------------------------------------
I've made an effort to make these portable to most Octave installations
as well as for use in most Matlab versions.  These work with no modification
in Matlab v5.2 &amp; v5.3.  If you want to use these in Matlab, however,
you'll do yourself a favor by renaming ode45.m and ode23.m to something else,
like ode45_octave.m and ode23_octave.m.  This is because Matlab already has
two integrators named ode45 and ode23.
Don't forget to change the function names at the top of ode45_octave.m and
ode23_octave.m as well.

Some effort has been made to create ode45 and ode23 with similar
argument structures to their Matlab counterparts.  Only the most basic function
calls to ode45 and ode23 match up in both Matlab's integrators and these.
Feel free to change them as you see fit.  You are welcome to mail me for useful
change suggestions.

----------------------------------------------------------------------
Basic differences between the integrators:

In general, the higher the integration order, the smaller the local truncation
error is at each time step.  Small local truncation errors result in larger
integration steps.  This is demonstrated by ode78 generating far fewer
steps than ode23 for solving the same problem over the same time interval with
the same error criterion.
The cost of the higher order integrators is the number of function evaluations
required at each step.  This results in longer execution times for each integration
step however, if the higher order integrator takes much larger steps, the
total computational cost is less.
The tradeoff between solving a problem with an integrator that takes fewer steps
versus using one that takes more time for each step will vary with each problem.
Factors such as numerical stiffness (or lack of) and the number of discontinuities
present in the ODE's will somteimes cause ode23 to be more effective than ode78,
&amp; vice-versa.

Integrator costs in right-hand-side (RHS) evaluations:
   - ode23.m    : requires 3 RHS function evaluations per step
   - ode45.m    : requires 6 RHS function evaluations per step (+1 for the Dormand-Prince pair)
   - ode78.m    : requires 13 RHS function evaluations per step

   - rk2fixed.m : requires 2 RHS function evaluations per step
   - rk4fixed.m : requires 4 RHS function evaluations per step
   - rk8fixed.m : requires 13 RHS function evaluations per step

If you are interested in experimenting with pendulum.m, try turning on
the 'trace' variable for screen output.  This increases execution time
but provides a way to monitor the problem during the integration.

----------------------------------------------------------------------
Two of the original files rk4fixed.m and ode78.m came from other
people.

rk4fixed.m was written by:
     Dr. Raul Longoria, Dept. of Mechanical Engineering,
     The Univ. of Texas at Austin and

ode78.m was originally written by:
     Dr. Howard Wilson &amp; 
     Daljeet Singh, Dept. Of Electrical Engineering,
     The University of Alabama

The other files were created from the structure of these 
along with coefficients from standard numerical methods books.
Dr. Longoria has given permission to redistribute rk4fixed.m.

ode78.m was originally found at:
ftp://ftp.mathworks.com/pub/contrib/v4/diffeq/ode78.m.
&amp; I am redistributing a modified version with Dr. Wilson's
permission.  

Numerous applications in ordinary and partial differential equations
can be found in Dr. Wilson's text:

     Howard Wilson and Louis Turcotte, 'Advanced Mathematics and 
     Mechanics Applications Using MATLAB', 2nd Ed, CRC Press, 1997

----------------------------------------------------------------------

Bugs or changes or comments should be directed to the email address
below.  The latest and greatest versions will generally _NOT_ be found
on the octave-source list, but rather will be maintained  at:
	http://marc.me.utexas.edu/tmp/octave_ode_solvers/
and hopefully sometime soon at:
	http://octave.sourceforge.net/


Marc Compere
CompereM at asme dot org
created : 06 October 1999
modified: 04 August 2002</PRE>

<BR>
---------------------------------------------------------
<BR>
<TT>ode23.m:</TT>
<BR>
-------
<PRE>function [tout,xout] = ode23(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% ode23.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% ode23.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% ode23 (v1.14) Integrates a system of ordinary differential equations using
% 2nd &amp; 3rd order Runge-Kutta formulas.  This particular 3rd-order method reduces
% to Simpson's 1/3 rule and uses the 3rd order estimate for xout.
%
% 3rd-order accurate RK methods have local and global errors of O(h^4) and O(h^3),
% respectively and yield exact results when the solution is a cubic.
%
% The order of the RK method is the order of the local *truncation* error, d,
% which is the principle error term in the portion of the Taylor series
% expansion that gets dropped, or intentionally truncated.  This is different
% from the local error which is the difference between the estimated solution
% and the actual, or true solution.  The local error is used in stepsize
% selection and may be approximated by the difference between two estimates of
% different order, l(h) = x_(O(h+1)) - x_(O(h)).  With this definition, the
% local error will be as large as the error in the lower order method.
% The local truncation error is within the group of terms that gets multipled
% by h when solving for a solution from the general RK method.  Therefore, the
% order-p solution created by the RK method will be roughly accurate to O(h^(p+1))
% since the local truncation error shows up in the solution as h*d, which is
% h times an O(h^(p)) term, or rather O(h^(p+1)).
% Summary:   For an order-p accurate RK method,
%            - the local truncation error is O(h^p)
%            - the local error used for stepsize adjustment and that
%              is actually realized in a solution is O(h^(p+1))
%
% This requires 3 function evaluations per integration step.
%
% Relevant discussion on step size choice can be found on pp.90,91 in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential Equations
% and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% The error estimate formula and slopes are from
% Numerical Methods for Engineers, 2nd Ed., Chapra &amp; Canale, McGraw-Hill, 1985
%
% Usage:
%         [tout, xout] = ode23(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)
%
% INPUT:
% FUN   - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where FUN = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - Initial value COLUMN-vector.
% ode_fcn_format - this specifies if the user-defined ode function is in
%         the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%         or:           xprime = fun(x,t)   (ode_fcn_format=1)
%         Matlab's solvers comply with ode_fcn_format=0 while
%         Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% tol   - The desired accuracy. (optional, default: tol = 1.e-6).
% trace - If nonzero, each step is printed. (optional, default: trace = 0).
% count - if nonzero, variable 'rhs_counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'rhs_counter' is incremented in here, not in the state-dot file
%         simply make 'rhs_counter' global in the file that calls ode23
% hmax  - limit the maximum stepsize to be less than or equal to hmax
%
% OUTPUT:
% tout  - Returned integration time points (column-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 27 June 2001

if nargin &lt; 8, hmax = (tspan(2) - tspan(1))/2.5; end
if nargin &lt; 7, count = 0; end
if nargin &lt; 6, trace = 0; end
if nargin &lt; 5, tol = 1.e-3; end
if nargin &lt; 4, ode_fcn_format = 0; end

pow = 1/4; % see p.91 in the Ascher &amp; Petzold reference for more infomation.

% The 2(3) coefficients:
 a(1,1)=0;
 a(2,1)=1/2;
 a(3,1)=-1; a(3,2)=2;
 % 2nd order b-coefficients
 b2(1)=0; b2(2)=1;
 % 5th order b-coefficients
 b3(1)=1/6; b3(2)=2/3; b3(3)=1/6;
 for i=1:3
  c(i)=sum(a(i,:));
 end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmin = (tfinal - t)/1e12;
h = (tfinal - t)/200; % initial guess at a step size
x = x0(:);            % this always creates a column vector, x
tout = t;             % first output time
xout = x.';           % first output solution

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0; end
end % if count

if trace
 clc, t, h, x
end

% The main loop
   while (t &lt; tfinal) &amp; (h &gt;= hmin)
      if t + h &gt; tfinal, h = tfinal - t; end

      % compute the slopes
      if (ode_fcn_format==0), % (default)
         k(:,1)=feval(FUN,t,x);
         k(:,2)=feval(FUN,t+c(2)*h,x+h*(a(2,1)*k(:,1)));
         k(:,3)=feval(FUN,t+c(3)*h,x+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)));
      else, % ode_fcn_format==1
         k(:,1)=feval(FUN,x,t);
         k(:,2)=feval(FUN,x+h*(a(2,1)*k(:,1)),t+c(2)*h);
         k(:,3)=feval(FUN,x+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)),t+c(3)*h);
      end % if (ode_fcn_format==1)

      % increment rhs_counter
      if count==1, rhs_counter = rhs_counter + 3; end

      % compute the 2nd order estimate
      x2=x + h*b2(2)*k(:,2);
      % compute the 3rd order estimate
      x3=x + h*(b3(1)*k(:,1) + b3(2)*k(:,2) + b3(3)*k(:,3));

      % estimate the local truncation error
      gamma1 = x3 - x2;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');
      tau = tol*max(norm(x,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if delta &lt;= tau
         t = t + h;
         x = x3;    % &lt;-- using the higher order estimate is called 'local extrapolation'
         tout = [tout; t];
         xout = [xout; x.'];
      end
      if trace
         home, t, h, x
      end

      % Update the step size
      if delta == 0.0
       delta = 1e-16;
      end
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end;

   if (t &lt; tfinal)
      disp('Step size grew too small.')
      t, h, x
   end</PRE>

<BR>
---------------------------------------------------------
<BR>
<TT>ode45.m:</TT>
<BR>
-------
<PRE>function [tout,xout,Nsteps_acc,Nsteps_rej] = ode45(FUN,tspan,x0,pair,ode_fcn_format,tol,trace,count,hmax,N_est_acc_steps)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% ode45.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% ode45.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% ode45 (v1.15) integrates a system of ordinary differential equations using
% 4th &amp; 5th order embedded formulas from Dormand &amp; Prince or Fehlberg.
%
% The Fehlberg 4(5) pair is established and works well, however, the
% Dormand-Prince 4(5) pair minimizes the local truncation error in the
% 5th-order estimate which is what is used to step forward (local extrapolation.)
% Generally it produces more accurate results and costs roughly the same
% computationally.  The Dormand-Prince pair is the default.
%
% This is a 4th-order accurate integrator therefore the local error normally
% expected is O(h^5).  However, because this particular implementation
% uses the 5th-order estimate for xout (i.e. local extrapolation) moving
% forward with the 5th-order estimate should yield local error of O(h^6).
%
% The order of the RK method is the order of the local *truncation* error, d,
% which is the principle error term in the portion of the Taylor series
% expansion that gets dropped, or intentionally truncated.  This is different
% from the local error which is the difference between the estimated solution
% and the actual, or true solution.  The local error is used in stepsize
% selection and may be approximated by the difference between two estimates of
% different order, l(h) = x_(O(h+1)) - x_(O(h)).  With this definition, the
% local error will be as large as the error in the lower order method.
% The local truncation error is within the group of terms that gets multipled
% by h when solving for a solution from the general RK method.  Therefore, the
% order-p solution created by the RK method will be roughly accurate to O(h^(p+1))
% since the local truncation error shows up in the solution as h*d, which is
% h times an O(h^(p)) term, or rather O(h^(p+1)).
% Summary:   For an order-p accurate RK method,
%            - the local truncation error is O(h^p)
%            - the local error used for stepsize adjustment and that
%              is actually realized in a solution is O(h^(p+1))
%
% This requires 6 function evaluations per integration step.
%
% Both the Dormand-Prince and Fehlberg 4(5) coefficients are from a tableu in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential Equations
% and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% The error estimate formula and slopes are from
% Numerical Methods for Engineers, 2nd Ed., Chappra &amp; Cannle, McGraw-Hill, 1985
%
% Usage:
%         [tout, xout] = ode45(FUN,tspan,x0,pair,ode_fcn_format,tol,trace,count,hmax)
%
% INPUTS:
%    FUN   - String containing name of user-supplied problem description.
%            Call: xprime = fun(t,x) where FUN = 'fun'.
%            t      - Time (scalar).
%            x      - Solution column-vector.
%            xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
%    tspan - [ tstart, tfinal ]
%    x0    - Initial value COLUMN-vector.
%    pair  - flag specifying which integrator coefficients to use:
%               0 --&gt; use Dormand-Prince 4(5) pair (default)
%               1 --&gt; use Fehlberg pair 4(5) pair
%    ode_fcn_format - this specifies if the user-defined ode function is in
%            the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%            or:           xprime = fun(x,t)   (ode_fcn_format=1)
%            Matlab's solvers comply with ode_fcn_format=0 while
%            Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
%    tol   - The desired accuracy. (optional, default: tol = 1.e-6).
%    trace - If nonzero, each step is printed. (optional, default: trace = 0).
%    count - if nonzero, variable 'rhs_counter' is initalized, made global
%            and counts the number of state-dot function evaluations
%            'rhs_counter' is incremented in here, not in the state-dot file
%            simply make 'rhs_counter' global in the file that calls ode45
%    hmax  - limit the maximum stepsize to be less than or equal to hmax
%    N_est_acc_steps - an estimate of how many accepted steps there may be;
%                      used to preallocate memory for the [tout,xout] solutions
%
% OUTPUTS:
%    tout  - array of discretized times points (an Nsteps_acc by 1 column-vector).
%    xout  - solution, one solution column-vector per tout-value (an Nsteps_acc by Nstates matrix)
%    Nsteps_acc - total number of accepted integration steps + 1
%    Nsteps_rej - total number of rejected integration steps
%
% The result can be displayed by: plot(tout, xout).
%
% The following relationship should hold after a completed run:
%    rhs_counter == (Nsteps_acc-1+Nsteps_rej)*6+1
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 03 July 2001

if nargin &lt; 10, N_est_acc_steps = (tspan(2)-tspan(1))*1e3; end
if nargin &lt;  9, hmax = (tspan(2) - tspan(1))/2.5; end
if nargin &lt;  8, count = 0; end
if nargin &lt;  7, trace = 0; end
if nargin &lt;  6, tol = 1.e-6; end
if nargin &lt;  5, ode_fcn_format = 0; end
if nargin &lt;  4, pair = 0; end

pow = 1/6; % see p.91 in the Ascher &amp; Petzold reference for more infomation.

if (pair==0),
   % Use the Dormand-Prince 4(5) coefficients:
    a_(1,1)=0;
    a_(2,1)=1/5;
    a_(3,1)=3/40; a_(3,2)=9/40;
    a_(4,1)=44/45; a_(4,2)=-56/15; a_(4,3)=32/9;
    a_(5,1)=19372/6561; a_(5,2)=-25360/2187; a_(5,3)=64448/6561; a_(5,4)=-212/729;
    a_(6,1)=9017/3168; a_(6,2)=-355/33; a_(6,3)=46732/5247; a_(6,4)=49/176; a_(6,5)=-5103/18656;
    a_(7,1)=35/384; a_(7,2)=0; a_(7,3)=500/1113; a_(7,4)=125/192; a_(7,5)=-2187/6784; a_(7,6)=11/84;
    % 4th order b-coefficients
    b4_(1,1)=5179/57600; b4_(2,1)=0; b4_(3,1)=7571/16695; b4_(4,1)=393/640; b4_(5,1)=-92097/339200; b4_(6,1)=187/2100; b4_(7,1)=1/40;
    % 5th order b-coefficients
    b5_(1,1)=35/384; b5_(2,1)=0; b5_(3,1)=500/1113; b5_(4,1)=125/192; b5_(5,1)=-2187/6784; b5_(6,1)=11/84; b5_(7,1)=0;
    for i=1:7,
     c_(i)=sum(a_(i,:));
    end
else, % pair==1 so use the Fehlberg 4(5) coefficients:
    a_(1,1)=0;
    a_(2,1)=1/4;
    a_(3,1)=3/32; a_(3,2)=9/32;
    a_(4,1)=1932/2197; a_(4,2)=-7200/2197; a_(4,3)=7296/2197;
    a_(5,1)=439/216; a_(5,2)=-8; a_(5,3)=3680/513; a_(5,4)=-845/4104;
    a_(6,1)=-8/27; a_(6,2)=2; a_(6,3)=-3544/2565; a_(6,4)=1859/4104; a_(6,5)=-11/40;
    % 4th order b-coefficients (guaranteed to be a column vector)
    b4_(1,1)=25/216; b4_(2,1)=0; b4_(3,1)=1408/2565; b4_(4,1)=2197/4104; b4_(5,1)=-1/5;
    % 5th order b-coefficients (also guaranteed to be a column vector)
    b5_(1,1)=16/135; b5_(2,1)=0; b5_(3,1)=6656/12825; b5_(4,1)=28561/56430; b5_(5,1)=-9/50; b5_(6,1)=2/55;
    for i=1:6,
     c_(i)=sum(a_(i,:));
    end
end % if (pair==0)


% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmin = (tfinal - t)/1e20;
h = (tfinal - t)/100; % initial step size guess
x = x0(:);            % ensure x is a column vector

Nstates = size(x,1);
tout = zeros(N_est_acc_steps,1);       % preallocating memory is not only faster but, for long
xout = zeros(N_est_acc_steps,Nstates); % simulations, prevents disappointing surprises later

Nsteps_rej = 0;
Nsteps_acc = 1;
tout(Nsteps_acc) = t;     % first output time
xout(Nsteps_acc,:) = x.'; % first output solution = IC's

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0; end
end % if count

if trace
 clc, t, h %, x
end

if (pair==0),

   k_ = x*zeros(1,7);  % k_ needs to be initialized as an Nx7 matrix where N=number of rows in x
                       % (just for speed so octave doesn't need to allocate more memory at each stage value)

   % Compute the first stage prior to the main loop.  This is part of the Dormand-Prince pair caveat.
   % Normally, during the main loop the last stage for x(k) is the first stage for computing x(k+1).
   % So, the very first integration step requires 7 function evaluations, then each subsequent step
   % 6 function evaluations because the first stage is simply assigned from the last step's last stage.
   % note: you can see this by the last element in c_ is 1.0, thus t+c_(7)*h = t+h, ergo, the next step.
   if (ode_fcn_format==0), % (default)
      k_(:,1)=feval(FUN,t,x); % first stage
   else, % ode_fcn_format==1
      k_(:,1)=feval(FUN,x,t);
   end % if (ode_fcn_format==1)

   % increment rhs_counter
   if (count==1), rhs_counter = rhs_counter + 1; end

   % The main loop using Dormand-Prince 4(5) pair
   while (t &lt; tfinal) &amp; (h &gt;= hmin),
      if t + h &gt; tfinal, h = tfinal - t; end

      % Compute the slopes by computing the k(:,j+1)'th column based on the previous k(:,1:j) columns
      % notes: k_ needs to end up as an Nxs, a_ is 7x6, which is s by (s-1),
      %        s is the number of intermediate RK stages on [t (t+h)] (Dormand-Prince has s=7 stages)
      if (ode_fcn_format==0), % (default)
         for j = 1:6,
            k_(:,j+1) = feval(FUN, t+c_(j+1)*h, x+h*k_(:,1:j)*a_(j+1,1:j)');
         end
      else, % ode_fcn_format==1
         for j = 1:6,
            k_(:,j+1) = feval(FUN, x+h*k_(:,1:j)*a_(j+1,1:j)', t+c_(j+1)*h);
         end
      end % if (ode_fcn_format==1)

      % increment rhs_counter
      if (count==1), rhs_counter = rhs_counter + 6; end

      % compute the 4th order estimate
      x4=x + h* (k_*b4_); % k_ is Nxs (or Nx7) and b4_ is a 7x1

      % compute the 5th order estimate
      x5=x + h*(k_*b5_); % k_ is Nxs (or Nx7) and b5_ is a 7x1

      % estimate the local truncation error
      gamma1 = x5 - x4;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');       % actual error
      tau = tol*max(norm(x,'inf'),1.0); % allowable error

      % Update the solution only if the error is acceptable
      if (delta&lt;=tau),

         t = t + h;
         x = x5;    % &lt;-- using the higher order estimate is called 'local extrapolation'
         Nsteps_acc = Nsteps_acc + 1;
         tout(Nsteps_acc) = t;
         xout(Nsteps_acc,:) = x.';

         % Assign the last stage for x(k) as the first stage for computing x(k+1).
         % This is part of the Dormand-Prince pair caveat.
         % k_(:,7) has already been computed, so use it instead of recomputing it
         % again as k_(:,1) during the next step.
         k_(:,1)=k_(:,7);

      else, % unacceptable integration step

         Nsteps_rej = Nsteps_rej + 1;

      end % if (delta&lt;=tau)

      if trace
         home, t, h, Nsteps_acc, Nsteps_rej
      end % if trace

      % Update the step size
      if (delta==0.0),
       delta = 1e-16;
      end % if (delta==0.0)
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end % while (t &lt; tfinal) &amp; (h &gt;= hmin)

else, % pair==1 --&gt; use RK-Fehlberg 4(5) pair

   k_ = x*zeros(1,6);  % k_ needs to be initialized as an Nx6 matrix where N=number of rows in x
                       % (just for speed so octave doesn't need to allocate more memory at each stage value)
 
   % The main loop using Dormand-Prince 4(5) pair
   while (t &lt; tfinal) &amp; (h &gt;= hmin),
      if t + h &gt; tfinal, h = tfinal - t; end

      % Compute the slopes by computing the k(:,j+1)'th column based on the previous k(:,1:j) columns
      % notes: k_ needs to end up as an Nx6, a_ is 6x5, which is s by (s-1),  (RK-Fehlberg has s=6 stages)
      %        s is the number of intermediate RK stages on [t (t+h)]
      if (ode_fcn_format==0), % (default)
         k_(:,1)=feval(FUN,t,x); % first stage
      else, % ode_fcn_format==1
         k_(:,1)=feval(FUN,x,t);
      end % if (ode_fcn_format==1)

      if (ode_fcn_format==0), % (default)
         for j = 1:5,
            k_(:,j+1) = feval(FUN, t+c_(j+1)*h, x+h*k_(:,1:j)*a_(j+1,1:j)');
         end
      else, % ode_fcn_format==1
         for j = 1:5,
            k_(:,j+1) = feval(FUN, x+h*k_(:,1:j)*a_(j+1,1:j)', t+c_(j+1)*h);
         end
      end % if (ode_fcn_format==1)

      % increment rhs_counter
      if (count==1), rhs_counter = rhs_counter + 6; end

      % compute the 4th order estimate
      x4=x + h* (k_(:,1:5)*b4_); % k_(:,1:5) is an Nx5 and b4_ is a 5x1

      % compute the 5th order estimate
      x5=x + h*(k_*b5_); % k_ is the same as k_(:,1:6) and is an Nx6 and b5_ is a 6x1

      % estimate the local truncation error
      gamma1 = x5 - x4;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');       % actual error
      tau = tol*max(norm(x,'inf'),1.0); % allowable error

      % Update the solution only if the error is acceptable
      if (delta&lt;=tau),

         t = t + h;
         x = x5;    % &lt;-- using the higher order estimate is called 'local extrapolation'
         Nsteps_acc = Nsteps_acc + 1;
         tout(Nsteps_acc) = t;
         xout(Nsteps_acc,:) = x.';

      else, % unacceptable integration step

         Nsteps_rej = Nsteps_rej + 1;

      end % if (delta&lt;=tau)

      if trace
         home, t, h, Nsteps_acc, Nsteps_rej
      end % if trace

      % Update the step size
      if (delta==0.0),
       delta = 1e-16;
      end % if (delta==0.0)
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end % while (t &lt; tfinal) &amp; (h &gt;= hmin)

end % if (pair==0),


% if necessary, trim outputs
if (Nsteps_acc&lt;N_est_acc_steps),
   tout = tout(1:Nsteps_acc);
   xout = xout(1:Nsteps_acc,:);
end % if (Nsteps_acc&lt;N_est_acc_steps)


if (t &lt; tfinal),
   disp('Step size grew too small.')
   t, h, Nsteps_acc, Nsteps_rej
end</PRE>

<BR>
---------------------------------------------------------
<BR>
<TT>ode78.m:</TT>
<BR>
-------
<PRE>function [tout,xout] = ode78(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% ode78.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% ode78.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% ode78 (v1.14) Integrates a system of ordinary differential equations using
% 7th order formulas.
%
% This is a 7th-order accurate integrator therefore the local error normally
% expected is O(h^8).  However, because this particular implementation
% uses the 8th-order estimate for xout (i.e. local extrapolation) moving
% forward with the 8th-order estimate will yield errors on the order of O(h^9).
%
% The order of the RK method is the order of the local *truncation* error, d,
% which is the principle error term in the portion of the Taylor series
% expansion that gets dropped, or intentionally truncated.  This is different
% from the local error which is the difference between the estimated solution
% and the actual, or true solution.  The local error is used in stepsize
% selection and may be approximated by the difference between two estimates of
% different order, l(h) = x_(O(h+1)) - x_(O(h)).  With this definition, the
% local error will be as large as the error in the lower order method.
% The local truncation error is within the group of terms that gets multipled
% by h when solving for a solution from the general RK method.  Therefore, the
% order-p solution created by the RK method will be roughly accurate to O(h^(p+1))
% since the local truncation error shows up in the solution as h*d, which is
% h times an O(h^(p)) term, or rather O(h^(p+1)).
% Summary:   For an order-p accurate RK method,
%            - the local truncation error is O(h^p)
%            - the local error used for stepsize adjustment and that
%              is actually realized in a solution is O(h^(p+1))
%
% This requires 13 function evaluations per integration step.
%
% Relevant discussion on step size choice can be found on pp.90,91 in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential Equations
% and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% More may be found in the original author's text containing numerous
% applications on ordinary and partial differential equations using Matlab:
%
%     Howard Wilson and Louis Turcotte, 'Advanced Mathematics and 
%     Mechanics Applications Using MATLAB', 2nd Ed, CRC Press, 1997
%
%
% [tout, xout] = ode78(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)
%
% INPUT:
% FUN   - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where FUN = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - Initial value COLUMN-vector.
% ode_fcn_format - this specifies if the user-defined ode function is in
%         the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%         or:           xprime = fun(x,t)   (ode_fcn_format=1)
%         Matlab's solvers comply with ode_fcn_format=0 while
%         Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% tol   - The desired accuracy. (optional, default: tol = 1.e-6).
% trace - If nonzero, each step is printed. (optional, default: trace = 0).
% count - if nonzero, variable 'rhs_counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'rhs_counter' is incremented in here, not in the state-dot file
%         simply make 'rhs_counter' global in the file that calls ode78
% hmax  - limit the maximum stepsize to be less than or equal to hmax
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).

%   Daljeet Singh &amp; Howard Wilson
%   Dept. Of Electrical Engg., The University of Alabama.
%   11-24-1988.
%
% modified by:
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 19 May 2001


% The Fehlberg coefficients:
alpha_ = [ 2./27., 1/9, 1/6, 5/12, 0.5, 5/6, 1/6, 2/3, 1/3, 1, 0, 1 ]';
beta_  = [ 2/27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/36, 1/12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/24, 0, 1/8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          5/12, 0, -25/16, 25/16, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          0.05, 0, 0, 0.25, 0.2, 0, 0, 0, 0, 0, 0, 0, 0 ;
          -25/108, 0, 0, 125/108, -65/27, 125/54, 0, 0, 0, 0, 0, 0, 0 ;
          31/300, 0, 0, 0, 61/225, -2/9, 13/900, 0, 0, 0, 0, 0, 0 ;
          2, 0, 0, -53/6, 704/45, -107/9, 67/90, 3, 0, 0, 0, 0, 0 ;
          -91/108, 0, 0, 23/108, -976/135, 311/54, -19/60, 17/6, -1/12, 0, 0, 0, 0 ;
          2383/4100, 0, 0, -341/164, 4496/1025, -301/82, 2133/4100, 45/82, 45/164, 18/41, 0, 0, 0 ;
          3/205, 0, 0, 0, 0, -6/41, -3/205, -3/41, 3/41, 6/41, 0, 0, 0 ;
          -1777/4100, 0, 0, -341/164, 4496/1025, -289/82, 2193/4100, 51/82, 33/164, 12/41, 0, 1, 0 ]';
chi_  = [ 0, 0, 0, 0, 0, 34/105, 9/35, 9/35, 9/280, 9/280, 0, 41/840, 41/840]';
psi_  = [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1 ]';

pow = 1/8; % see p.91 in the Ascher &amp; Petzold reference for more infomation.


if nargin &lt; 8, hmax = (tspan(2) - tspan(1))/2.5; end
if nargin &lt; 7, count = 0; end
if nargin &lt; 6, trace = 0; end
if nargin &lt; 5, tol = 1.e-6; end
if nargin &lt; 4, ode_fcn_format = 0; end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
% the following step parameters are used in ODE45
% hmax = (tfinal - t)/5;
% hmin = (tfinal - t)/20000;
% h = (tfinal - t)/100;
% The following parameters were taken because the integrator has
% higher order than ODE45. This choice is somewhat subjective.
%hmin = (tfinal - t)/10000;
hmin = (tfinal - t)/1e20;
h = (tfinal - t)/50;
x = x0(:);          % the '(:)' ensures x is initialized as a column vector
f = x*zeros(1,13);  % f needs to be an Nx13 matrix where N=number of rows in x
tout = t;
xout = x.';
tau = tol * max(norm(x,'inf'), 1);

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

if trace
%  clc, t, h, x
%   clc, t, x
   clc, t
end

% The main loop
   while (t &lt; tfinal) &amp; (h &gt;= hmin)
      if t + h &gt; tfinal, h = tfinal - t; end

      % Compute the slopes
      if (ode_fcn_format==0), % (default)
         f(:,1) = feval(FUN,t,x);
         for j = 1: 12,
            f(:,j+1) = feval(FUN, t+alpha_(j)*h, x+h*f*beta_(:,j));
         end
      else, % ode_fcn_format==1
         f(:,1) = feval(FUN,x,t);
         for j = 1: 12,
            f(:,j+1) = feval(FUN, x+h*f*beta_(:,j), t+alpha_(j)*h);
         end
      end % if (ode_fcn_format==1)


      % increment rhs_counter
      if count==1, rhs_counter = rhs_counter + 13; end

      % Truncation error term
      gamma1 = h*41/840*f*psi_;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');
      tau = tol*max(norm(x,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if delta &lt;= tau
         t = t + h;
         x = x + h*f*chi_;  % this integrator uses local extrapolation
         tout = [tout; t];
         xout = [xout; x.'];
      end
      if trace
%        home, t, h, x
%        home, t, x
         home, t, h
      end

      % Update the step size
      if delta == 0.0
       delta = 1e-16;
      end
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end;

   if (t &lt; tfinal)
      disp('SINGULARITY LIKELY.')
      t
   end</PRE>

<BR>
---------------------------------------------------------
<BR>
<TT>penddot.m:</TT>
<BR>
---------
<PRE>function xdot=penddot(t,x)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% penddot.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% penddot.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% This is an example derivative function file that works
% with Octave or Matlab.
% The equations represent motion of a simple pendulum with damping.
%
% The plots created by pendulum.m show the angular position and velocity
% trajectories created by each different integrator.
% Position is the trace that reaches steady state at 0(rad)
% because of the gravity term, -m*g*l/2*sin(x(1)).
% Velocity reaches a steady state of zero because of the
% damping term, -b*x(2).
%
% Use ode45 to integrate these ODE's
% like this:
%    [t,x] = ode45('penddot',tspan,IC);
%
% x is the state column vector and meant to be used only within this m-file
% This function is meant to return the derivatives of the state variable
% given the state vector and time.
%
% Structure:	xdot = [ x1dot, x2dot, ..., xNdot ]'
%
% eg.  ml^2*thetadd + b*thetad + m*g*l*sin(theta) = 0
%
%	x(1) = theta
%	x(2) = thetad 	( = x(1)dot )
%
% Convention: the lowest order states are first columnwise
%
% Marc Compere
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 23 October 2001

global m g l b counter index

xdot=[ x(2) , 1/(1/3*m*l^2)*(-b*x(2)-m*g*l/2*sin(x(1)))]';

% remember to return a column vector

end</PRE>

<BR>
---------------------------------------------------------
<BR>
<TT>pendulum.m:</TT>
<BR>
----------
<PRE>% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% pendulum.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% pendulum.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% This integrates a set of ordinary differential equations (ODE) using 6 different
% ODE solvers.  The equations represent the dynamics of a simple pendulum.
%
% The integrators ode78.m, ode45.m, ode23.m, rk8fixed.m, rk4fixed.m, and rk2fixed.m
% all produce column vector output similar to Matlab.
%
% All integrators work in octave 2.1.32 and Matlab 5.3 with no modification.
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 23 October 2001

clear

% allow global access to the parameters m, g, l, &amp; b:
global m g l b rhs_counter

m=1;    % (kg)
g=9.81; % (m/s^2)
l=1.;   % (m)
b=0.7;   % ((N-s)/m))

% integrator setup:
trace = 0;   % this is a (1/0) flag that puts output to the screen or not
count = 1;   % this is a (1/0) flag that causes rk4fixed to increment 'rhs_counter' or not
rhs_counter = 0; % 'rhs_counter' is the number of right-hand-side function evaluations, z_dot=f(z)
ode_fcn_format = 0; % 0 chooses Matlab-format right-hand-side function definitions
                    % xdot=f(t,x), or Octave's lsode format, xdot=f(x,t)
t0=0;
tfinal = 5; % (seconds)
tspan = [t0,tfinal];
hmax = 0.1;

% Initial Conditions: theta(t=0)=30(deg) &amp; initially at rest
IC = [ 30*pi/180,0]'; % (rad), (rad/s)
sps = 100;            % sps -&gt; step per second
Nsteps=(tfinal-t0)*sps;    % this creates sps number of integration steps per second
tolerance = 1e-3;

% polite housekeeping
this_user_uses=page_screen_output;
page_screen_output=0;

% for development:
tol=tolerance; pair=0; x0=IC; FUN='penddot';

% Solve the ODE specified in penddot.m using each of the 6 m-file integrators.

t_begin_calcs=cputime;

   disp('Integrating using rk2fixed...')
   [t1,zrk2fixed] = rk2fixed('penddot',tspan,IC,Nsteps,ode_fcn_format,trace,count);     % fixed step integration
   rk2_counter = rhs_counter;
   rhs_counter=0;

t(1)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using rk4fixed...')
   [t2,zrk4fixed] = rk4fixed('penddot',tspan,IC,Nsteps,ode_fcn_format,trace,count);     % fixed step integration
   rk4_counter = rhs_counter;
   rhs_counter=0;

t(2)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using rk8fixed...')
   [t3,zrk8fixed] = rk8fixed('penddot',tspan,IC,Nsteps,ode_fcn_format,trace,count);     % fixed step integration
   rk8_counter = rhs_counter;
   rhs_counter=0;

t(3)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode23...')
   [t4,zode23] = ode23('penddot',tspan,IC,ode_fcn_format,tolerance,trace,count,hmax); % rk23 variable step integration
   ode23_counter = rhs_counter;
   rhs_counter=0;

t(4)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode45 with the Dormand-Prince 4(5) pair...')
   pair=0;
   [t5,zode45dp,Nsteps_acc_ode45dp,Nsteps_rej_ode45dp] = ode45('penddot',tspan,IC,pair,ode_fcn_format,tolerance,trace,count,hmax); % rk45 variable step integration
   ode45dp_counter = rhs_counter;
   rhs_counter=0;

t(5)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode45 with the Fehlberg 4(5) pair...')
   pair=1;
   [t6,zode45rkf,Nsteps_acc_ode45rkf,Nsteps_rej_ode45rkf] = ode45('penddot',tspan,IC,pair,ode_fcn_format,tolerance,trace,count,hmax); % rk45 variable step integration
   ode45rkf_counter = rhs_counter;
   rhs_counter=0;

t(6)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode78...')
   [t7,zode78] = ode78('penddot',tspan,IC,ode_fcn_format,tolerance,trace,count,hmax); % rk78 variable step integration
   ode78_counter = rhs_counter;
   rhs_counter=0;

t(7)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   %disp('Integrating using sdirk...')
   % Note: If you want to use sdirk you have to compile sdirk.oct, then change the function definition in
   %       penddot.m to zdot=penddot(z,t) and uncomment the lines below.
   %       Then if you want the other integrators to still work, set ode_fcn_format=1 above.
   %skip_step=10;
   %h_initial = 1e-3;
   %[t8,zsdirk] = sdirk('penddot',IC,tspan,tolerance,0.1*tolerance,h_initial,trace,skip_step); % 4th order stiff integration
   %sdirk_counter = rhs_counter;
   %rhs_counter=0;

t(8)=cputime-t_begin_calcs;

rhs_fcn_evaluation_summary=[
strcat('rk2_counter      = ',num2str(rk2_counter));
strcat('rk4_counter      = ',num2str(rk4_counter));
strcat('rk8_counter      = ',num2str(rk8_counter));
strcat('ode23_counter    = ',num2str(ode23_counter));
strcat('ode45dp_counter  = ',num2str(ode45dp_counter));
strcat('ode45rkf_counter = ',num2str(ode45rkf_counter));
strcat('ode78_counter    = ',num2str(ode78_counter)) ]
%strcat('sdirk_counter    = ',num2str(sdirk_counter)) ]


disp('Elapsed times for each solver to integrate the state equations for a simple pendulum:')
t


% notes: The Dormand-Prince pair in ode45 produces a time output vector
%        of 18x1.  (size(t5)=18x1)
%        The number of function evaluations is 103 which you can compute
%        from (18-1)*6+1=103.  (18-1) because there were really only 17 new
%        steps computed.  The first step is just the initial conditions.
%        Multiply (18-1) by 6 because each step during the normal main loop
%        requires 6 function evaluations to compute.  Add 1 because the very
%        first step requires 1 function evaluation to start the main loop.


% plot those puppies
figure(1)
clg
if ishold~=1, hold, end
title('Pendulum Position &amp; Velocity')
ylabel('Theta &amp; Theta_dot (rad) &amp; (rad/s)')
xlabel('Time (s)')
plot(t1,zrk2fixed)
plot(t2,zrk4fixed)
plot(t3,zrk8fixed)
plot(t4,zode23)
plot(t5,zode45dp)
plot(t6,zode45rkf)
plot(t7,zode78)
%plot(t8,zsdirk)
if ishold==1, hold, end


% These plots show the angular position and velocity
% trajectories created by each different integrator.
% Position is the trace that reaches steady state at -pi/2.
% Velocity reaches a steady state of zero.


% Return setting(s) to what they were before
page_screen_output=this_user_uses;</PRE>

<BR>

<BR>
---------------------------------------------------------
<BR>
<TT>rk2fixed.m:</TT>
<BR>
----------
<PRE>function [tout,xout] = rk2fixed(FUN,tspan,x0,Nsteps,ode_fcn_format,trace,count)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% rk2fixed.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% rk2fixed.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% rk2fixed (v1.14) integrates a system of ordinary differential equations using a
% 2nd order Runge-Kutta formula called Ralston's method.
% This choice of 2nd order coefficients provides a minimum bound on truncation error.
% For more, see Ralston &amp; Rabinowitz (1978) or 
% Numerical Methods for Engineers, 2nd Ed., Chappra &amp; Cannle, McGraw-Hill, 1985
%
% 2nd-order accurate RK methods have a local error estimate of O(h^3).
%
% rk2fixed() requires 2 function evaluations per integration step.
%
% Usage:
%         [tout, xout] = rk2fixed(FUN, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% FUN    - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where FUN = 'fun'.
%          t      - Time (scalar).
%          x      - Solution column-vector.
%          xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan  - [ tstart, tfinal ]
% x0     - Initial value COLUMN-vector.
% Nsteps - number of steps used to span [ tstart, tfinal ]
% ode_fcn_format - this specifies if the user-defined ode function is in
%          the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%          or:           xprime = fun(x,t)   (ode_fcn_format=1)
%          Matlab's solvers comply with ode_fcn_format=0 while
%          Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% trace  - If nonzero, each step is printed. (optional, default: trace = 0).
% count  - if nonzero, variable 'rhs_counter' is initalized, made global
%          and counts the number of state-dot function evaluations
%          'rhs_counter' is incremented in here, not in the state-dot file
%          simply make 'rhs_counter' global in the file that calls rk2fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 19 May 2001

if nargin &lt; 7, count = 0; end
if nargin &lt; 6, trace = 0; end
if nargin &lt; 5, Nsteps = 400/(tspan(2)-tspan(1)); end % &lt;-- 400 is a guess for a default,
                                                %  try verifying the solution with rk4fixed
if nargin &lt; 4, ode_fcn_format = 0; end

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/Nsteps;
xout(1,:) = x0';           
tout(1)   = t;
x = x0(:);

if trace
 clc, t, h, x
end

% The main loop
h = (tspan(2)-tspan(1))/Nsteps;

for i=1:Nsteps,
     if (ode_fcn_format==0),
      k1 = feval(FUN,t,x);
      k2 = feval(FUN,t+3/4*h,x+3/4*h*k1);
     else,
      k1 = feval(FUN,x,t);
      k2 = feval(FUN,x+3/4*h*k1,t+3/4*h);
     end % if (ode_fcn_format==0)

     % increment rhs_counter
     if count==1,
      rhs_counter = rhs_counter + 2;
     end % if

     t = t + h;
     x = (x+h*(1/3*k1+2/3*k2));
     tout = [tout; t];
     xout = [xout; x.'];

     if trace,
      home, t, h, x
     end

end</PRE>

<BR>
---------------------------------------------------------
<BR>
<TT>rk4fixed.m:</TT>
<BR>
----------
<PRE>function [tout,xout] = rk4fixed(FUN,tspan,x0,Nsteps,ode_fcn_format,trace,count)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% rk4fixed.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% rk4fixed.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% rk4fixed (v1.14) is a 4th order Runge-Kutta numerical integration routine.
% It requires 4 function evaluations per step.
%
% 4th-order accurate RK methods have a local error estimate of O(h^5).
%
%
% Usage:
%         [tout, xout] = rk4fixed(FUN, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% FUN    - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where FUN = 'fun'.
%          t      - Time or independent variable (scalar).
%          x      - Solution column-vector.
%          xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan  - [ tstart, tfinal ]
% x0     - Initial value COLUMN-vector.
% Nsteps - number of steps used to span [ tstart, tfinal ]
% ode_fcn_format - this specifies if the user-defined ode function is in
%          the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%          or:           xprime = fun(x,t)   (ode_fcn_format=1)
%          Matlab's solvers comply with ode_fcn_format=0 while
%          Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% trace  - If nonzero, each step is printed. (optional, default: trace = 0).
% count  - if nonzero, variable 'rhs_counter' is initalized, made global
%          and counts the number of state-dot function evaluations
%          'rhs_counter' is incremented in here, not in the state-dot file
%          simply make 'rhs_counter' global in the file that calls rk4fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 19 May 2001

if nargin &lt; 7, count = 0; end
if nargin &lt; 6, trace = 0; end
if nargin &lt; 5, Nsteps = 100/(tspan(2)-tspan(1)); end % &lt;-- 100 is a guess for a default,
                                                %  try verifying the solution with rk8fixed
if nargin &lt; 4, ode_fcn_format = 0; end

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/Nsteps;
xout(1,:) = x0';
tout(1) = t;
x = x0(:);
halfh = 0.5*h;

if trace
 clc, t, h, x
end

for i=1:Nsteps,
     if (ode_fcn_format==0),
      RK1 = feval(FUN,t,x);
      thalf = t+halfh;
      xtemp = x+halfh*RK1;
      RK2 = feval(FUN,thalf,xtemp);
      xtemp = x+halfh*RK2;
      RK3 = feval(FUN,thalf,xtemp);
      tfull = t+h;
      xtemp = x+h*RK3;
      RK4 = feval(FUN,tfull,xtemp);
     else,
      RK1 = feval(FUN,x,t);
      thalf = t+halfh;
      xtemp = x+halfh*RK1;
      RK2 = feval(FUN,xtemp,thalf);
      xtemp = x+halfh*RK2;
      RK3 = feval(FUN,xtemp,thalf);
      tfull = t+h;
      xtemp = x+h*RK3;
      RK4 = feval(FUN,xtemp,tfull);
     end % if (ode_fcn_format==0)

     % increment rhs_counter
     if count==1,
      rhs_counter = rhs_counter + 4;
     end % if

     t = t + h;
     x = (x+h/6*(RK1+2.0*(RK2+RK3)+RK4));
     tout = [tout; t];
     xout = [xout; x.'];

     if trace,
      home, t, h, x
     end

end</PRE>

<BR>

<BR>
---------------------------------------------------------
<BR>
<TT>rk8fixed.m:</TT>
<BR>
----------
<PRE>function [tout,xout] = rk8fixed(FUN,tspan,x0,Nsteps,ode_fcn_format,trace,count)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% rk8fixed.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% rk8fixed.m is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% rk8fixed (v1.14) is an 8th order Runge-Kutta numerical integration routine.
% It requires 13 function evaluations per step.  This is not the most
% efficient 8th order implementation.  It was just the easiest to put
% together as a variant from ode78.m.
%
% 8th-order accurate RK methods have a local error estimate of O(h^9).
%
%
% Usage:
%         [tout, xout] = rk8fixed(FUN, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% FUN    - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where FUN = 'fun'.
%          t      - Time or independent variable (scalar).
%          x      - Solution column-vector.
%          xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan  - [ tstart, tfinal ]
% x0     - Initial value COLUMN-vector.
% Nsteps - number of steps used to span [ tstart, tfinal ]
% ode_fcn_format - this specifies if the user-defined ode function is in
%          the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%          or:           xprime = fun(x,t)   (ode_fcn_format=1)
%          Matlab's solvers comply with ode_fcn_format=0 while
%          Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% trace  - If nonzero, each step is printed. (optional, default: trace = 0).
% count  - if nonzero, variable 'rhs_counter' is initalized, made global
%          and counts the number of state-dot function evaluations
%          'rhs_counter' is incremented in here, not in the state-dot file
%          simply make 'rhs_counter' global in the file that calls rk4fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM at asme dot org
% created : 06 October 1999
% modified: 19 May 2001

if nargin &lt; 7, count = 0; end
if nargin &lt; 6, trace = 0; end
if nargin &lt; 5, Nsteps = 50/(tspan(2)-tspan(1)); end % &lt;-- 50 is a guess for a default,
                                                %  try verifying the solution with ode78
if nargin &lt; 4, ode_fcn_format = 0; end

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

alpha_ = [ 2./27., 1/9, 1/6, 5/12, 0.5, 5/6, 1/6, 2/3, 1/3, 1, 0, 1 ]';
beta_  = [ 2/27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/36, 1/12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/24, 0, 1/8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          5/12, 0, -25/16, 25/16, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          0.05, 0, 0, 0.25, 0.2, 0, 0, 0, 0, 0, 0, 0, 0 ;
          -25/108, 0, 0, 125/108, -65/27, 125/54, 0, 0, 0, 0, 0, 0, 0 ;
          31/300, 0, 0, 0, 61/225, -2/9, 13/900, 0, 0, 0, 0, 0, 0 ;
          2, 0, 0, -53/6, 704/45, -107/9, 67/90, 3, 0, 0, 0, 0, 0 ;
          -91/108, 0, 0, 23/108, -976/135, 311/54, -19/60, 17/6, -1/12, 0, 0, 0, 0 ;
          2383/4100, 0, 0, -341/164, 4496/1025, -301/82, 2133/4100, 45/82, 45/164, 18/41, 0, 0, 0 ;
          3/205, 0, 0, 0, 0, -6/41, -3/205, -3/41, 3/41, 6/41, 0, 0, 0 ;
          -1777/4100, 0, 0, -341/164, 4496/1025, -289/82, 2193/4100, 51/82, 33/164, 12/41, 0, 1, 0 ]';
chi_   = [ 0, 0, 0, 0, 0, 34/105, 9/35, 9/35, 9/280, 9/280, 0, 41/840, 41/840]';

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/Nsteps;
xout(1,:) = x0';
tout(1) = t;
x = x0(:);
f = x*zeros(1,13);

if trace
 clc, t, h, x
end

for i=1:Nsteps,

     % Compute the slopes
     if (ode_fcn_format==0),
      f(:,1) = feval(FUN,t,x);
      for j = 1:12
         f(:,j+1) = feval(FUN, t+alpha_(j)*h, x+h*f*beta_(:,j));
      end
     else,
      f(:,1) = feval(FUN,x,t);
      for j = 1:12
         f(:,j+1) = feval(FUN, x+h*f*beta_(:,j), t+alpha_(j)*h);
      end
     end % if (ode_fcn_format==0)

     % increment rhs_counter
     if count==1,
      rhs_counter = rhs_counter + 13;
     end % if

     t = t + h;
     x = x + h*f*chi_;
     tout = [tout; t];
     xout = [xout; x.'];

     if trace,
      home, t, h, x
     end

end</PRE>

<BR>
---------------------------------------------------------
<BR>

<BR>

<BR>

<BR>

<TABLE CELLSPACING="0" CELLPADDING="0" WIDTH="100%">
<TR>
<TD>
<PRE>-- 
_______________________________________________________
 Marc Compere, Ph.D., The University of Texas at Austin
 CompereM at asme dot org, (512) 826-0729,   &lt;&gt;&lt;
 http://marc.me.utexas.edu/
_______________________________________________________</PRE>
</TD>
</TR>
</TABLE>

</BODY>
</HTML>

--=-+D8QRHJhVFcPnLLnycDs--