From octave-sources-request at bevo dot che dot wisc dot edu  Mon May 15 12:37:17 2000
Subject: ode solvers v1.07
From: Marc Compere <compere at mail dot utexas dot edu>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Mon, 15 May 2000 12:35:21 -0500

This is a multi-part message in MIME format.
--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Octave users,

	This message's attachments are files that comprise the contents of
http://marc.me.utexas.edu/tmp/octave_ode_solvers/ode_v1.07.tar.gz.  They
are attached in the following order:

      ode23.m
      ode45.m
      ode78.m
      rk2fixed.m
      rk4fixed.m
      rk8fixed.m
      pendulum.m
      penddot.m
      readme.txt

Best regards,
Marc Compere

-- 
_________________________________________________
Marc Compere, The University of Texas at Austin
e-mail: compere at mail dot utexas dot edu
www   : http://nerdlab.me.utexas.edu/~compere
office: ETC 4.134a, (512) 471-7347
_________________________________________________
--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="ode23.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="ode23.m"

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

% Copyright (C) 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.07) Integrates a system of ordinary differential equations using
% 2nd & 3rd order Runge-Kutta formulas.  The particular 3rd order method is
% Simpson's 1/3 rule.
% This particular implementation uses the 3rd order estimate for xout, although
% the truncation error is of order(h^2), therefore this method, overall is
% a 2nd-order method.
% This requires 3 function evaluations per integration step.
%
% The error estimate formula and slopes are from
% Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle, McGraw-Hill, 1985
%
% Usage:
%         [tout, xout] = ode23(F, tspan, x0, ode_fcn_format, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where F = '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
%
% 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
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 15 May 2000

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/8;

% 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;
hmax = (tfinal - t)/2.5;
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),
       k(:,1)=feval(F,t,x);
       k(:,2)=feval(F,t+c(2)*h,x+h*(a(2,1)*k(:,1)));
       k(:,3)=feval(F,t+c(3)*h,x+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)));
      else,
       k(:,1)=feval(F,x,t);
       k(:,2)=feval(F,x+h*(a(2,1)*k(:,1)),t+c(2)*h);
       k(:,3)=feval(F,x+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)),t+c(3)*h);
      end % if (ode_fcn_format==0)

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

      % 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


--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="ode45.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="ode45.m"

function [tout, xout] = ode45(F,tspan,x0,ode_fcn_format,tol,trace,count)

% Copyright (C) 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.07) integrates a system of ordinary differential equations using
% 4th & 5th order embedded Runge-Kutta-Fehlberg formulas.
% This requires 6 function evaluations per integration step.
% This particular implementation uses the 5th order estimate for xout, although
% the truncation error is of order(h^4), therefore this method, overall is
% a 4th-order method.
%
% The 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(F, tspan, x0, ode_fcn_format, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where F = '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 ode45
%
% 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
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 15 May 2000

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

pow = 1/8;

% 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
 b4_(1)=25/216; b4_(2)=0; b4_(3)=1408/2565; b4_(4)=2197/4104; b4_(5)=-1/5;
 % 5th order b-coefficients
 b5_(1)=16/135; b5_(2)=0; b5_(3)=6656/12825; b5_(4)=28561/56430; b5_(5)=-9/50; b5_(6)=2/55;
 for i=1:6
  c_(i)=sum(a_(i,:));
 end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmax = (tfinal - t)/2.5;
hmin = (tfinal - t)/1e9;
h = (tfinal - t)/100; % 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),
       k_(:,1)=feval(F,t,x);
       k_(:,2)=feval(F,t+c_(2)*h,x+h*(a_(2,1)*k_(:,1)));
       k_(:,3)=feval(F,t+c_(3)*h,x+h*(a_(3,1)*k_(:,1)+a_(3,2)*k_(:,2)));
       k_(:,4)=feval(F,t+c_(4)*h,x+h*(a_(4,1)*k_(:,1)+a_(4,2)*k_(:,2)+a_(4,3)*k_(:,3)));
       k_(:,5)=feval(F,t+c_(5)*h,x+h*(a_(5,1)*k_(:,1)+a_(5,2)*k_(:,2)+a_(5,3)*k_(:,3)+a_(5,4)*k_(:,4)));
       k_(:,6)=feval(F,t+c_(6)*h,x+h*(a_(6,1)*k_(:,1)+a_(6,2)*k_(:,2)+a_(6,3)*k_(:,3)+a_(6,4)*k_(:,4)+a_(6,5)*k_(:,5)));
      else,
       k_(:,1)=feval(F,x,t);
       k_(:,2)=feval(F,x+h*(a_(2,1)*k_(:,1)),t+c_(2)*h);
       k_(:,3)=feval(F,x+h*(a_(3,1)*k_(:,1)+a_(3,2)*k_(:,2)),t+c_(3)*h);
       k_(:,4)=feval(F,x+h*(a_(4,1)*k_(:,1)+a_(4,2)*k_(:,2)+a_(4,3)*k_(:,3)),t+c_(4)*h);
       k_(:,5)=feval(F,x+h*(a_(5,1)*k_(:,1)+a_(5,2)*k_(:,2)+a_(5,3)*k_(:,3)+a_(5,4)*k_(:,4)),t+c_(5)*h);
       k_(:,6)=feval(F,x+h*(a_(6,1)*k_(:,1)+a_(6,2)*k_(:,2)+a_(6,3)*k_(:,3)+a_(6,4)*k_(:,4)+a_(6,5)*k_(:,5)),t+c_(6)*h);
      end % if (ode_fcn_format==0)

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

      % compute the 4th order estimate
      x4=x + h*(b4_(1)*k_(:,1) + b4_(3)*k_(:,3) + b4_(4)*k_(:,4) + b4_(5)*k_(:,5));
      % compute the 5th order estimate
      x5=x + h*(b5_(1)*k_(:,1) + b5_(3)*k_(:,3) + b5_(4)*k_(:,4) + b5_(5)*k_(:,5) + b5_(6)*k_(:,6));

      % 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'
         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


--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="ode78.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="ode78.m"

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

% Copyright (C) 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.07) Integrates a system of ordinary differential equations using
% 7th order formulas.
% This particular implementation uses the 8th order estimate for xout, which
% is called 'local extrapolation'.  The truncation error, gamma1, is of order(h^8).
% Therefore this method, overall is a 7th-order method.
% This requires 13 function evaluations per integration step.
%
% 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(F, tspan, x0, ode_fcn_format, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where F = '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
%
% 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
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 15 May 2000


% The Fehlberg coefficients:
alpha_ = [ 2./27. 1/9 1/6 5/12 .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  ]
[ .05   0  0  .25  .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;

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.
hmax = (tfinal - t)/2.5;
%hmin = (tfinal - t)/10000;
hmin = (tfinal - t)/1000000000;
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 cols 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
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),
       f(:,1) = feval(F,t,x);
       for j = 1: 12
          f(:,j+1) = feval(F, t+alpha_(j)*h, x+h*f*beta_(:,j));
       end
      else,
       f(:,1) = feval(F,x,t);
       for j = 1: 12
          f(:,j+1) = feval(F, 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

      % 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
      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


--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="rk2fixed.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="rk2fixed.m"

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

% Copyright (C) 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.07) 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
%
% rk2fixed() requires 2 function evaluations per integration step.
%
% Usage:
%         [tout, xout] = rk2fixed(F, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% F      - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where F = '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
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 15 May 2000

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(F,t,x);
      k2 = feval(F,t+3/4*h,x+3/4*h*k1);
     else,
      k1 = feval(F,x,t);
      k2 = feval(F,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


--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="rk4fixed.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="rk4fixed.m"

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

% Copyright (C) 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.07) is a 4th order Runge-Kutta numerical integration routine.
% It requires 4 function evaluations per step.
%
% Usage:
%         [tout, xout] = rk4fixed(F, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% F      - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where F = '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
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 15 May 2000

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;

for i=1:Nsteps,
     if (ode_fcn_format==0),
      RK1 = feval(F,t,x);
      thalf = t+halfh;
      xtemp = x+halfh*RK1;
      RK2 = feval(F,thalf,xtemp);
      xtemp = x+halfh*RK2;
      RK3 = feval(F,thalf,xtemp);
      tfull = t+h;
      xtemp = x+h*RK3;
      RK4 = feval(F,tfull,xtemp);
     else,
      RK1 = feval(F,x,t);
      thalf = t+halfh;
      xtemp = x+halfh*RK1;
      RK2 = feval(F,xtemp,thalf);
      xtemp = x+halfh*RK2;
      RK3 = feval(F,xtemp,thalf);
      tfull = t+h;
      xtemp = x+h*RK3;
      RK4 = feval(F,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


--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="rk8fixed.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="rk8fixed.m"

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

% Copyright (C) 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.07) 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.
%
% Usage:
%         [tout, xout] = rk8fixed(F, tspan, x0, Nsteps, ode_fcn_format, trace, count)
%
% INPUT:
% F      - String containing name of user-supplied problem derivatives.
%          Call: xprime = fun(t,x) where F = '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
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 15 May 2000

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 .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  ]
[ .05   0  0  .25  .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);

for i=1:Nsteps,

     % Compute the slopes
     if (ode_fcn_format==0),
      f(:,1) = feval(F,t,x);
      for j = 1:12
         f(:,j+1) = feval(F, t+alpha_(j)*h, x+h*f*beta_(:,j));
      end
     else,
      f(:,1) = feval(F,x,t);
      for j = 1:12
         f(:,j+1) = feval(F, 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


--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="pendulum.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="pendulum.m"

% Copyright (C) 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.24 and Matlab 5.3 with no modification.
%
% Marc Compere
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 15 May 2000

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=2.;   % ((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];
% 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;

% 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); % rk45 variable step integration
   ode23_counter = rhs_counter
   rhs_counter=0;

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

   disp('Integrating using ode45...')
   [t5,zode45] = ode45('penddot',tspan,IC,ode_fcn_format,tolerance,trace,count); % rk45 variable step integration
   ode45_counter = rhs_counter
   rhs_counter=0;

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

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

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

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

t(7)=cputime-t_begin_calcs;



disp('Elapsed times for each solver to integrate a pendulum:')
t



% plot that baby
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,zode45)
plot(t6,zode78)
%plot(t7,zsdirk4)
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.



--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="penddot.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="penddot.m"

function zdot=penddot(t,z)

% Copyright (C) 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 -pi/2
% because of the gravity term, -m*g*l/2*cos(z(1)).
% Velocity reaches a steady state of zero because of the
% damping term, -b*z(2).
%
% Use ode45 to integrate these ODE's
% like this:
%    [t,z] = ode45('penddot',tspan,IC);
%
% z 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:	zdot = [ z1dot, z2dot, ... zNdot ]'
%
% eg.  ml^2*thetadd + b*thetad + m*g*l*sin(theta) = 0
%
%	z(1) = theta
%	z(2) = thetad 	( = z(1)dot )
%
% Convention: the lowest order states are first columnwise
%
% Marc Compere
% compere at mail dot utexas dot edu
% created : 06 October 1999
% modified: 15 May 2000

global m g l b counter index

zdot=[ z(2) , 1/(1/3*m*l^2)*(-b*z(2)-m*g*l/2*cos(z(1)))]';

% remember to return a column vector

end


--------------FA3BCD0FC824BBCC2CD89D5F
Content-Type: text/plain; charset=us-ascii;
 name="readme.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="readme.txt"

This is the readme.txt file for Octave m-file ordinary differential
equation (ODE) solvers, version 1.07.  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.
All are explicit RK formulas that work well with nonstiff problems.

----------------------------------------------------------------------
The archive octave_ode_solvers_v1.07.tar.gz contains 6 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.14 or better),
from a unix shell:

(1) unzip and untar the archive:
	gunzip ode_v1.07.tar.gz
	tar xvf ode_v1.07.tar
(2) change directories into the newly created directory and execute Octave:
	cd ode_v1.07
	octave
(3) run the sample pemdulum script from within Octave.  This will sequentially
    execute all 6 m-file integrators:
	pendulum

----------------------------------------------------------------------
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.
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 the degree of numerical stiffness or the number of discontinuities
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
   - 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 speed
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 origianlly 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/

Marc Compere
compere at mail dot utexas dot edu
created : 06 October 1999
modified: 15 May 2000


--------------FA3BCD0FC824BBCC2CD89D5F--



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------