--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,13 @@
+2002-07-10  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* ODE.h: Move integrate and do_integrate method declarations and
+	definitions here.
+	* base-de.h: From here.
+
+	* ODES.h, ODES.cc, ODESFunc.h, ODESSA.h, ODESSA.cc: New files.
+	* Makefile.in: Add them to the appropriate lists.
+	(LIBOCTAVE_CXX_SOURCES): 
+
 2002-07-02  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* NLEqn.cc (NLEqn::error_message): New function.

--- a/liboctave/Makefile.in
+++ b/liboctave/Makefile.in
@@ -46,7 +46,8 @@
 
 INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.h DASPK.h DASSL.h FEGrid.h \
 	LinConst.h LP.h LPsolve.h LSODE.h NLConst.h NLEqn.h NLFunc.h \
-	NLP.h ODE.h ODEFunc.h Objective.h QP.h Quad.h Range.h base-de.h \
+	NLP.h ODE.h ODEFunc.h ODES.h ODESFunc.h ODESSA.h \
+	Objective.h QP.h Quad.h Range.h base-de.h \
 	base-min.h byte-swap.h cmd-edit.h cmd-hist.h data-conv.h \
 	dir-ops.h file-ops.h file-stat.h getopt.h glob-match.h \
 	idx-vector.h lo-ieee.h lo-mappers.h lo-specfun.h lo-sysdep.h \
@@ -84,9 +85,10 @@
 	vx-rv-cs.cc vx-s-ccv.cc vx-s-crv.cc \
 	vx-rv-crv.cc vx-cv-ccv.cc vx-crv-rv.cc vx-ccv-cv.cc
 
-LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc DAE.cc DASPK.cc \
-	DASSL.cc FEGrid.cc LinConst.cc LPsolve.cc LSODE.cc \
-	NLEqn.cc Quad.cc Range.cc data-conv.cc dir-ops.cc \
+LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc DAE.cc \
+	DASPK.cc DASSL.cc FEGrid.cc LinConst.cc LPsolve.cc \
+	LSODE.cc NLEqn.cc ODES.cc ODESSA.cc Quad.cc \
+	Range.cc data-conv.cc dir-ops.cc \
 	file-ops.cc file-stat.cc glob-match.cc idx-vector.cc \
 	lo-ieee.cc lo-mappers.cc lo-specfun.cc lo-sysdep.cc \
 	lo-utils.cc mach-info.cc oct-alloc.cc oct-env.cc \

--- a/liboctave/ODE.h
+++ b/liboctave/ODE.h
@@ -51,6 +51,64 @@
     }
 
   ~ODE (void) { }
+
+  // Derived classes must provide functions to actually do the
+  // integration.
+
+  // Return the vector of states at output time t.
+  virtual ColumnVector do_integrate (double tt) = 0;
+
+  // Return a matrix of states at each output time specified by t.
+  // The rows of the result matrix should each correspond to a new
+  // output time.
+  virtual Matrix do_integrate (const ColumnVector& tt) = 0;
+
+  virtual Matrix do_integrate (const ColumnVector& tt,
+			       const ColumnVector& ttcrit) = 0;
+
+  // Lots of ways to call the single function and optionally set and
+  // get additional information.
+
+  // Integrate to t from current point.
+  virtual ColumnVector integrate (double tt)
+    { return do_integrate (tt); }
+
+  // Set new x0, t0 and integrate to t.
+  virtual ColumnVector integrate (const ColumnVector& x0, double t0, double tt)
+    {
+      initialize (x0, t0);
+      return do_integrate (tt);
+    }
+
+  // Integrate from current point and return output at all points
+  // specified by t.
+  virtual Matrix integrate (const ColumnVector& tt)
+    { return do_integrate (tt); }
+
+  // Set new x0, t0 and integrate to return output at all points
+  // specified by t.
+  virtual Matrix integrate (const ColumnVector& x0, double t0,
+			    const ColumnVector& tt)
+    {
+      initialize (x0, t0);
+      return do_integrate (tt);
+    }
+
+  // Integrate from current point and return output at all points
+  // specified by t.
+  virtual Matrix integrate (const ColumnVector& tt,
+			    const ColumnVector& ttcrit)
+    { return do_integrate (tt, ttcrit); }
+
+  // Set new x0, t0 and integrate to return output at all points
+  // specified by t.
+  virtual Matrix integrate (const ColumnVector& x0, double t0,
+			    const ColumnVector& tt,
+			    const ColumnVector& ttcrit)
+    {
+      initialize (x0, t0);
+      return do_integrate (tt, ttcrit);
+    }
 };
 
 #endif

new file mode 100644
--- /dev/null
+++ b/liboctave/ODES.cc
@@ -0,0 +1,54 @@
+/*
+
+Copyright (C) 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave 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.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "ODES.h"
+#include "lo-error.h"
+
+void
+ODES::initialize (const ColumnVector& xx, double t)
+{
+  base_diff_eqn::initialize (xx, t);
+  xdot = ColumnVector (xx.length (), 0.0);
+}
+
+void
+ODES::initialize (const ColumnVector& xx, double t,
+		  const ColumnVector& xtheta)
+{
+  base_diff_eqn::initialize (xx, t);
+  xdot = ColumnVector (xx.length (), 0.0);
+  theta = xtheta;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/ODES.h
@@ -0,0 +1,89 @@
+/*
+
+Copyright (C) 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave 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.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if !defined (octave_ODES_h)
+#define octave_ODES_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#include "ODESFunc.h"
+#include "base-de.h"
+
+class
+ODES : public base_diff_eqn, public ODESFunc
+{
+public:
+
+  ODES (void)
+    : base_diff_eqn (), ODESFunc (), theta () { }
+
+  ODES (const ColumnVector& x, double t, ODESFunc& f)
+    : base_diff_eqn (x, t), ODESFunc (f), xdot (x.length (), 0.0), theta () { }
+
+  ODES (const ColumnVector& x, const ColumnVector& xtheta, double t,
+	ODESFunc& f)
+    : base_diff_eqn (x, t), ODESFunc (f), xdot (x.length (), 0.0),
+      theta (xtheta) { }
+
+  ODES (const ODES& a)
+    : base_diff_eqn (a), ODESFunc (a), theta (a.theta) { }
+
+  ODES& operator = (const ODES& a)
+    {
+      if (this != &a)
+	{
+	  base_diff_eqn::operator = (a);
+	  ODESFunc::operator = (a);
+
+	  xdot = a.xdot;
+	  theta = a.theta;
+	}
+      return *this;
+    }
+
+  ~ODES (void) { }
+
+  ColumnVector parameter_vector (void) { return theta; }
+
+  void initialize (const ColumnVector& x, double t);
+
+  void initialize (const ColumnVector& x, double t,
+		   const ColumnVector& theta);
+
+protected:
+
+  // State vector time derivatives.
+  ColumnVector xdot;
+
+  // Parameter vector.
+  ColumnVector theta;
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/ODESFunc.h
@@ -0,0 +1,117 @@
+/*
+
+Copyright (C) 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave 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.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if !defined (octave_ODESFunc_h)
+#define octave_ODESFunc_h 1
+
+#include "dMatrix.h"
+
+class
+ODESFunc
+{
+public:
+
+  struct DAEJac
+    {
+      Matrix *dfdxdot;
+      Matrix *dfdx;
+    };
+
+  typedef ColumnVector (*ODES_fsub) (double, const ColumnVector& x,
+				     const ColumnVector& theta); 
+
+  typedef ColumnVector (*ODES_bsub) (double, const ColumnVector& x,
+				     const ColumnVector& theta, int column);
+
+  typedef Matrix (*ODES_jsub) (double, const ColumnVector& x,
+			       const ColumnVector& theta);
+
+  ODESFunc (void)
+    : fsub (0), bsub (0), jsub (0) { }
+
+  ODESFunc (ODES_fsub f)
+    : fsub (f), bsub (0), jsub (0) { }
+
+  ODESFunc (ODES_fsub f, ODES_bsub b)
+    : fsub (f), bsub (b), jsub (0) { }
+
+  ODESFunc (ODES_fsub f, ODES_bsub b, ODES_jsub j)
+    : fsub (f), bsub (b), jsub (j) { }
+
+  ODESFunc (const ODESFunc& a)
+    : fsub (a.fsub), bsub (a.bsub), jsub (a.jsub) { }
+
+  ODESFunc& operator = (const ODESFunc& a)
+    {
+      if (this != &a)
+	{
+	  fsub = a.fsub;
+	  bsub = a.bsub;
+	  jsub = a.jsub;
+	}
+      return *this;
+    }
+
+  ~ODESFunc (void) { }
+
+  ODES_fsub fsub_function (void) const { return fsub; }
+
+  ODESFunc& set_fsub_function (ODES_fsub f)
+    {
+      fsub = f;
+      return *this;
+    }
+
+  ODES_bsub bsub_function (void) const { return bsub; }
+
+  ODESFunc& set_bsub_function (ODES_bsub b)
+    {
+      bsub = b;
+      return *this;
+    }
+
+  ODES_jsub jsub_function (void) const { return jsub; }
+
+  ODESFunc& set_jsub_function (ODES_jsub j)
+    {
+      jsub = j;
+      return *this;
+    }
+
+protected:
+
+  ODES_fsub fsub;
+  ODES_bsub bsub;
+  ODES_jsub jsub;
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
+
+
+
+

new file mode 100644
--- /dev/null
+++ b/liboctave/ODESSA.cc
@@ -0,0 +1,722 @@
+/*
+
+Copyright (C) 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave 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.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <iostream.h>
+
+#include <cstdlib>
+#include <cfloat>
+#include <cmath>
+
+// For instantiating the Array<Matrix> object.
+#include "Array.h"
+#include "Array.cc"
+
+#include "ODESSA.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "ov-fcn.h"
+#include "pager.h"
+#include "parse.h"
+#include "unwind-prot.h"
+#include "utils.h"
+#include "variables.h"
+
+#ifndef F77_FUNC
+#define F77_FUNC(x, X) F77_FCN (x, X)
+#endif
+void
+ODESSA_options::set_integration_method (const std::string& val)
+{
+  if (val == "stiff" || val == "bdf")
+    x_integration_method = "stiff";
+  else if (val == "non-stiff" || val == "adams")
+    x_integration_method = "non-stiff";
+  else
+    (*current_liboctave_error_handler)
+      ("odessa_options: method must be \"stiff\", \"bdf\", \"non-stiff\", or \"adams\"");
+}
+
+typedef int (*odessa_fcn_ptr) (int*, const double&, double*,
+                              double*, double*);
+
+typedef int (*odessa_jac_ptr) (int*, const double&, double*,
+                              double*, const int&, const int&,
+                              double*, const int&);
+
+typedef int (*odessa_dfdp_ptr) (int*, const double&, double*,
+                              double*, double*, const int&);
+
+
+extern "C"
+int F77_FUNC (odessa, ODESSA) (odessa_fcn_ptr, odessa_dfdp_ptr, int*,
+                               double*, double*, double&, double&,
+                               int&, double&, const double*, int&, 
+                               int&, int*, double*, int&, int*, int&,
+                               odessa_jac_ptr, int&);
+
+template class Array<Matrix>;
+
+static ODESFunc::ODES_fsub user_fsub;
+static ODESFunc::ODES_bsub user_bsub;
+static ODESFunc::ODES_jsub user_jsub;
+
+
+static int
+odessa_f (int* neq, const double& t, double *state,
+	  double *par, double *fval)
+{
+  int n = neq[0];
+  int n_par = neq[1];
+
+  // Load the state and parameter arrays as Octave objects
+
+  ColumnVector tmp_state (n);
+  for (int i = 0; i < n; i++)
+    tmp_state(i) = state[i];
+
+  ColumnVector tmp_param (n_par);
+  for (int i = 0; i < n_par; i++)
+    tmp_param(i) = par[i];
+
+  ColumnVector tmp_fval = user_fsub (t, tmp_state, tmp_param);
+
+  // Return the function value as a C array object
+
+  for (int i = 0; i < n; i++)
+    {
+      fval[i] = tmp_fval(i);
+    }
+
+  return 0;
+}
+
+static int
+odessa_j (int* neq, const double& t, double *state,
+	  double *par, const int& ml, const int& mu, double *pd, const
+	  int& nrowpd)
+{
+  int n = neq[0];
+  int n_par = neq[1];
+
+  // Load the state and parameter arrays as Octave objects
+  ColumnVector tmp_state (n);
+  for (int i = 0; i < n; i++)
+    tmp_state(i) = state[i];
+
+  ColumnVector tmp_param (n_par);
+  for (int i = 0; i < n_par; i++)
+    tmp_param(i) = par[i];
+
+  Matrix tmp_fval = user_jsub (t, tmp_state, tmp_param);
+
+  for (int j = 0; j < n; j++)
+    for (int i = 0; i < nrowpd; i++)
+      pd[nrowpd*j+i] = tmp_fval(i,j);
+
+  return 0;
+}
+
+static int
+odessa_b (int* neq, const double& t, double *state,
+	  double *par, double *dfdp, const int& jpar)
+
+{
+  int n = neq[0];
+  int n_par = neq[1];
+
+  // Load the state and parameter arrays as Octave objects
+  ColumnVector tmp_state (n);
+  for (int i = 0; i < n; i++)
+    tmp_state(i) = state[i];
+
+  ColumnVector tmp_param (n_par);
+  for (int i = 0; i < n_par; i++)
+    tmp_param(i) = par[i];
+
+  ColumnVector tmp_fval = user_bsub (t, tmp_state, tmp_param, jpar);
+
+  for (int i = 0; i < n; i++)
+    dfdp[i] = tmp_fval(i);
+
+  return 0;
+}
+
+ODESSA::ODESSA (void) : ODES (), ODESSA_options ()
+{
+  neq.resize(2);
+  n = size ();
+
+  stop_time_set = false;
+  stop_time = 0.0;
+
+  integration_error = false;
+  restart = true;
+  iopt.resize(4);
+
+  istate = 1;
+  itask = 1;
+  iopt(0) = 0;
+  isopt = 0;
+  iopt(1) = isopt;
+  npar = 0;
+  neq(0) = n;
+  neq(1) = npar;
+
+  sanity_checked = false;
+}
+
+ODESSA::ODESSA (const ColumnVector& state, double time, ODESFunc& f)
+  : ODES (state, time, f), ODESSA_options ()
+{
+  neq.resize(2);
+  n = size ();
+
+  stop_time_set = false;
+  stop_time = 0.0;
+ 
+  ODESFunc::operator = (f);
+
+  integration_error = false;
+  restart = true;
+
+  iopt.resize(4);
+  istate = 1;
+  itask = 1;
+  iopt(0) = 0;
+  isopt = 0;
+  iopt(1) = isopt;
+
+  sanity_checked = false;
+
+  npar = 0;
+  neq(0) = n;
+  neq(1) = npar;
+
+  y.resize (n, 1, 0.0);
+}
+
+ODESSA::ODESSA (const ColumnVector& state, const ColumnVector& theta,
+		const Matrix& sensitivity_guess, double time, ODESFunc& f)
+  : ODES (state, theta, time, f)
+{
+  initialized = false;
+  restart = false;
+
+  neq.resize(2);
+  n = state.length();
+  npar = theta.length();
+
+  neq(0) = n;
+  neq(1) = npar;
+
+  stop_time_set = false;
+  stop_time = 0.0;
+
+  ODESFunc::operator = (f);
+
+  sx0 = sensitivity_guess;
+  par.resize (npar);
+
+  for (int i = 0; i < npar; i++)
+    {
+      par(i) = theta(i);
+    }
+
+  sanity_checked = false;
+
+  npar = theta.length ();
+
+  iopt.resize(4);
+  istate = 1;
+  itask = 1;
+  iopt(0) = 0;
+  isopt = 1;
+  iopt(1) = isopt;
+
+  y.resize (n, npar+1, 0.0);
+}
+
+void
+ODESSA::force_restart (void)
+{
+  restart = true;
+}
+
+void
+ODESSA::set_stop_time (double time)
+{
+  stop_time_set = true;
+  stop_time = time;
+}
+
+void
+ODESSA::clear_stop_time (void)
+{
+  stop_time_set = 0;
+}
+
+void
+ODESSA::integrate (double tout)
+{
+  ODESSA_result retval;
+  if (! initialized)
+    {
+      
+      for (int i = 0; i < n; i++)
+	y(i,0) = x(i);
+      
+      if (npar > 0)
+	{
+	  for (int j = 0; j < npar; j++)
+	    for (int i = 0; i < n; i++)
+	      y(i,j+1) = sx0(i,j);
+	}
+      
+      integration_error = false;
+      
+      user_fsub = ODESFunc::fsub_function ();
+      user_bsub = ODESFunc::bsub_function ();
+      user_jsub = ODESFunc::jsub_function ();
+      
+      int idf;
+
+      if (user_bsub)
+	idf = 1;
+      else
+	idf = 0;
+      
+      iopt(2) = idf;
+      
+      
+      if (restart)
+	{
+	  restart = false;
+	  istate = 1;
+	}
+      
+      if (integration_method () == "stiff")
+	{
+	  if (user_jsub)
+	    method_flag = 21;
+	  else
+	    method_flag = 22;
+	  if (isopt)
+	    {
+	      liw = 21 + n + npar;
+	      lrw = 22 + 8*(npar+1)*n + n*n + n;
+	    }
+	  else
+	    {
+	      liw = 20 + n;
+	      lrw = 22 + 9*n + n*n;
+	    }
+	}
+      else
+	{
+	  if (isopt)
+	    {
+	      if (user_jsub)
+		method_flag = 11;
+	      else
+		method_flag = 12;
+	      liw = 21 + n + npar;
+	      lrw = 22 + 15*(npar+1)*n + n*n + n;
+	    }
+	  else
+	    {
+	      method_flag = 10;
+	      liw = 20 + n;
+	      lrw = 22 + 16 * n;
+	    }
+	}
+      
+      if (iwork.length () != liw)
+	{
+	  iwork.resize (liw);
+	  
+	  for (int i = 4; i < 10; i++)
+	    iwork.elem (i) = 0;
+	}
+      
+      if (rwork.length () != lrw)
+	{
+	  rwork.resize (lrw);
+	  
+	  for (int i = 4; i < 10; i++)
+	    rwork.elem (i) = 0.0;
+	}
+      
+      initialized = true;
+    }
+  integration_error = false;
+
+  // NOTE: this won't work if LSODE passes copies of the state vector.
+  //       In that case we have to create a temporary vector object
+  //       and copy.
+
+
+  if (! sanity_checked)
+    {
+      ColumnVector fval = user_fsub (t, x, theta);
+      
+      if (fval.length () != x.length ())
+	{
+	  (*current_liboctave_error_handler)
+	    ("odessa: inconsistent sizes for state and residual vectors");
+	  
+	  integration_error = true;
+	  return;
+	}
+
+      sanity_checked = true;
+    }
+
+  if (stop_time_set)
+    {
+      itask = 4;
+      rwork.elem (0) = stop_time;
+      iopt(0) = 1;
+    }
+  else
+    {
+      itask = 1;
+    }
+
+  double rel_tol = relative_tolerance ();
+
+  Array<double> abs_tol = absolute_tolerance ();  //note; this should
+  //  be a matrix, not a vector
+
+  int abs_tol_len = abs_tol.length ();
+
+  int itol;
+
+  if (abs_tol_len == 1)
+    itol = 1;
+  else if (abs_tol_len == n)
+    itol = 2;
+  else
+    {
+      (*current_liboctave_error_handler)
+        ("lsode: inconsistent sizes for state and absolute tolerance vectors");
+
+      integration_error = 1;
+      return;
+    }
+
+  if (initial_step_size () >= 0.0)
+    {
+      rwork.elem (4) = initial_step_size ();
+      iopt(0) = 1;
+    }
+
+  if (maximum_step_size () >= 0.0)
+    {
+      rwork.elem (5) = maximum_step_size ();
+      iopt(0) = 1;
+    }
+
+  if (minimum_step_size () >= 0.0)
+    {
+      rwork.elem (6) = minimum_step_size ();
+      iopt(0) = 1;
+    }
+
+  if (step_limit () > 0)
+    {
+      iwork.elem (5) = step_limit ();
+      iopt(0) = 1;
+    }
+
+      py = y.fortran_vec ();
+      piwork = iwork.fortran_vec ();
+      prwork = rwork.fortran_vec ();
+      ppar = par.fortran_vec ();
+      piopt = iopt.fortran_vec ();
+      pneq = neq.fortran_vec ();
+
+  const double *pabs_tol = abs_tol.fortran_vec ();
+
+   F77_XFCN (odessa, ODESSA, (odessa_f,odessa_b,pneq,py,ppar,t, tout,
+			      itol,rel_tol,pabs_tol,itask,istate,
+			      piopt,prwork,lrw,piwork,liw,odessa_j,
+			      method_flag));
+
+  if (f77_exception_encountered)
+    {
+      integration_error = true;
+      (*current_liboctave_error_handler) ("unrecoverable error in odessa");
+    }
+  else
+    {
+      switch (istate)
+        {
+        case -13: // Return requested in user-supplied function.
+        case -6:  // error weight became zero during problem. (solution
+                  // component i vanished, and atol or atol(i) = 0.)
+        case -5:  // repeated convergence failures (perhaps bad jacobian
+                  // supplied or wrong choice of mf or tolerances).
+        case -4:  // repeated error test failures (check all inputs).
+        case -3:  // illegal input detected (see printed message).
+        case -2:  // excess accuracy requested (tolerances too small).
+        case -1:  // excess work done on this call (perhaps wrong mf).
+          integration_error = 1;
+          break;
+	  
+	case 1:  // tout is same as t
+        case 2:  // lsode was successful
+          t = tout;
+          break;
+          
+        default: // Error?
+          (*current_liboctave_error_handler)
+            ("unrecognized value of istate returned from lsode");
+          break;
+        }
+    }
+}
+
+std::string
+ODESSA::error_message (void) const
+{
+  std::string retval;
+
+  switch (istate)
+    {
+    case -13:
+      retval = "return requested in user-supplied function";
+      break;
+
+    case -6:
+      retval = "Error weight became zero during problem.\
+  (solution component i vanished, and atol or atol(i) == 0)";
+      break;
+
+    case -5:
+      retval = "repeated convergence failures (perhaps bad jacobian\
+ supplied or wrong choice of integration method or tolerances)";
+      break;
+
+    case -4:
+      retval = "repeated error test failures (check all inputs)";
+      break;
+
+    case -3:
+      retval = "invalid input detected (see printed message)";
+      break;
+
+    case -2:
+      retval = "excess accuracy requested (tolerances too small)";
+      break;
+
+    case -1:
+      retval = "excess work on this call (perhaps wrong integration method)";
+      break;
+
+    case 1:
+      retval = "prior to initial call";
+      break;
+
+    case 2:
+      retval = "successful exit";
+      break;
+          
+    case 3:
+      retval = "prior to continuation call with modified parameters";
+      break;
+          
+    default:
+      retval = "unknown error state";
+      break;
+    }
+
+  return retval;
+}
+
+
+ODESSA_result
+ODESSA::integrate (const ColumnVector& tout)
+{
+  ODESSA_result retval;
+
+  Matrix x_out;
+
+  Array<Matrix> x_s_out;
+
+  int n_out = tout.capacity ();
+
+  if (n_out > 0 && n > 0)
+    {
+      x_out.resize (n_out, n);
+
+      x_s_out.resize (npar, Matrix (n_out, n, 0.0));
+
+      for (int j = 0; j < n_out; j++)
+	{
+	  integrate (tout(j));
+
+	  if (integration_error)
+	    {
+	      retval = ODESSA_result (x_out, x_s_out);
+	      return retval;
+	    }
+
+	  for (int i = 0; i < n; i++)
+	    {
+	      x_out(j,i) = y(i,0);
+
+	      for (int k = 0; k < npar; k++)
+		{
+		  x_s_out(k)(j,i) = y(i,k+1);
+		}
+	    }
+	}
+    }
+
+  retval = ODESSA_result (x_out, x_s_out);
+
+  return retval;
+}
+
+ODESSA_result
+ODESSA::integrate (const ColumnVector& tout, const ColumnVector& tcrit) 
+{
+  ODESSA_result retval;
+
+  Matrix x_out;
+
+  Array<Matrix> x_s_out;
+
+  int n_out = tout.capacity ();
+
+  if (n_out > 0 && n > 0)
+    {
+      x_out.resize (n_out, n);
+
+      x_s_out.resize (npar, Matrix (n_out, n, 0.0));
+
+      int n_crit = tcrit.capacity ();
+
+      if (n_crit > 0)
+	{
+	  int i_crit = 0;
+	  int i_out = 0;
+	  double next_crit = tcrit(0);
+	  double next_out;
+	  while (i_out < n_out)
+	    {
+	      bool do_restart = false;
+
+	      next_out = tout(i_out);
+	      if (i_crit < n_crit)
+		next_crit = tcrit(i_crit);
+
+	      int save_output;
+	      double t_out;
+
+	      if (next_crit == next_out)
+		{
+		  set_stop_time (next_crit);
+		  t_out = next_out;
+		  save_output = true;
+		  i_out++;
+		  i_crit++;
+		  do_restart = true;
+		}
+	      else if (next_crit < next_out)
+		{
+		  if (i_crit < n_crit)
+		    {
+		      set_stop_time (next_crit);
+		      t_out = next_crit;
+		      save_output = false;
+		      i_crit++;
+		      do_restart = true;
+		    }
+		  else
+		    {
+		      clear_stop_time ();
+		      t_out = next_out;
+		      save_output = true;
+		      i_out++;
+		    }
+		}
+	      else
+		{
+		  set_stop_time (next_crit);
+		  t_out = next_out;
+		  save_output = true;
+		  i_out++;
+		}
+	      integrate (t_out);
+	      if (integration_error)
+		{
+		  retval = ODESSA_result (x_out, x_s_out);
+		  return retval;
+		}
+	      if (save_output)
+		{
+		  for (int i = 0; i < n; i++)
+		    {
+		      x_out(i_out-1,i) = y(i,0);
+
+		      for (int k = 0; k < npar; k++)
+			{
+			  x_s_out(k)(i_out-1,i) = y(i,k+1);
+			}
+		    }
+		}
+
+	      if (do_restart)
+		force_restart ();
+	    }
+
+	  retval = ODESSA_result (x_out, x_s_out);
+	}
+      else
+	{
+	  retval = integrate (tout);
+
+	  if (integration_error)
+	    return retval;
+	}
+    }
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/ODESSA.h
@@ -0,0 +1,265 @@
+/*
+
+Copyright (C) 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave 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.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if !defined (octave_ODESSA_h)
+#define octave_ODESSA_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#include <cfloat>
+#include <cmath>
+
+#include "ODES.h"
+
+class
+ODESSA_options
+{
+public:
+
+  ODESSA_options (void) { init (); }
+
+  ODESSA_options (const ODESSA_options& opt) { copy (opt); }
+
+  ODESSA_options& operator = (const ODESSA_options& opt)
+    {
+      if (this != &opt)
+	copy (opt);
+
+      return *this;
+    }
+
+  ~ODESSA_options (void) { }
+
+  void init (void)
+    {
+      double sqrt_eps = ::sqrt (DBL_EPSILON);
+      x_absolute_tolerance.resize (1);
+      x_absolute_tolerance(0) = sqrt_eps;
+      x_initial_step_size = -1.0;
+      x_integration_method = "stiff";
+      x_maximum_step_size = -1.0;
+      x_minimum_step_size = 0.0;
+      x_relative_tolerance = sqrt_eps;
+
+      // This is consistent with earlier versions of Octave, and is
+      // much larger than the default of 500 specified in the LSODE
+      // sources.
+      x_step_limit = 100000;
+    }
+
+  void copy (const ODESSA_options& opt)
+    {
+      x_absolute_tolerance = opt.x_absolute_tolerance;
+      x_initial_step_size = opt.x_initial_step_size;
+      x_integration_method = opt.x_integration_method;
+      x_maximum_step_size = opt.x_maximum_step_size;
+      x_minimum_step_size = opt.x_minimum_step_size;
+      x_relative_tolerance = opt.x_relative_tolerance;
+      x_step_limit = opt.x_step_limit;
+    }
+
+  void set_default_options (void) { init (); }
+
+  void set_absolute_tolerance (double val)
+    {
+      x_absolute_tolerance.resize (1);
+      x_absolute_tolerance(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+    }
+
+  void set_absolute_tolerance (const Array<double>& val)
+    { x_absolute_tolerance = val; }
+
+  void set_initial_step_size (double val)
+    { x_initial_step_size = (val >= 0.0) ? val : -1.0; }
+
+  void set_integration_method (const std::string& val);
+
+
+  void set_maximum_step_size (double val)
+    { x_maximum_step_size = (val >= 0.0) ? val : -1.0; }
+
+  void set_minimum_step_size (double val)
+    { x_minimum_step_size = (val >= 0.0) ? val : 0.0; }
+
+  void set_relative_tolerance (double val)
+    { x_relative_tolerance = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); }
+
+  void set_step_limit (int val)
+    { x_step_limit = val; }
+
+  Array<double> absolute_tolerance (void) const
+    { return x_absolute_tolerance; }
+
+  double initial_step_size (void) const
+    { return x_initial_step_size; }
+
+  std::string integration_method (void) const
+    { return x_integration_method; }
+
+  double maximum_step_size (void) const
+    { return x_maximum_step_size; }
+
+  double minimum_step_size (void) const
+    { return x_minimum_step_size; }
+
+  double relative_tolerance (void) const
+    {  return x_relative_tolerance; }
+
+  int step_limit (void) const
+    { return x_step_limit; }
+
+private:
+
+  Array<double> x_absolute_tolerance;
+  double x_initial_step_size;
+  std::string x_integration_method;
+  double x_maximum_step_size;
+  double x_minimum_step_size;
+  double x_relative_tolerance;
+
+  int x_step_limit;
+};
+
+
+class
+ODESSA_result
+{
+public:
+
+  ODESSA_result (void) { }
+
+  ODESSA_result (const Matrix& xx, 
+		 const Array<Matrix>& xx_s)
+
+    : x (xx), x_s (xx_s) { }
+
+  ODESSA_result (const ODESSA_result& r)
+    : x (r.x), x_s (r.x_s) { }
+
+  ODESSA_result& operator = (const ODESSA_result& r)
+    {
+      if (this != &r)
+	{
+	  x = r.x;
+	  x_s = r.x_s;
+	}
+      return *this;
+    }
+
+  ~ODESSA_result (void) { }
+
+  Matrix state (void) const { return x; }
+  Array<Matrix> state_sensitivity (void) const { return x_s; }
+
+private:
+
+  Matrix x;
+  Array<Matrix> x_s;
+};
+
+class
+ODESSA : public ODES, public ODESSA_options
+{
+public:
+
+  ODESSA (void);
+
+  ODESSA (const ColumnVector& x, double time, ODESFunc& f);
+
+  ODESSA (const ColumnVector& x, const ColumnVector& theta,
+	  const Matrix& sensitivity_guess, double time, ODESFunc& f);
+
+  ~ODESSA (void) { }
+
+  void force_restart (void);
+
+  void set_stop_time (double t);
+  void clear_stop_time (void);
+
+  ODESSA_result integrate (const ColumnVector& tout);
+
+  ODESSA_result integrate (const ColumnVector& tout,
+			   const ColumnVector& tcrit); 
+
+  int integration_state (void) const { return istate; }
+
+  bool integration_ok (void) const { return ! integration_error; }
+
+  std::string error_message (void) const;
+
+private:
+
+  bool initialized;
+
+  bool sanity_checked;
+
+  bool stop_time_set;
+  double stop_time;
+
+  bool restart;
+
+  bool integration_error;
+
+  int liw;  
+  int lrw;
+  int method_flag;
+  Array<int> iwork;
+  Array<double> rwork;
+  int istate;
+  int itask;
+  Array<int> iopt;
+  int isopt;
+
+  Array<int> neq;
+
+  int n;
+  int npar;
+
+  // Hey, check out this crap: ZZZZ
+  Array<double> par;
+
+  Matrix sx0;
+
+  Matrix y;
+
+
+  double *py;
+  double *ppar;
+  int *piwork;
+  int *piopt;
+  int *pneq;
+  double *prwork;
+
+  void init_work_size (int);
+
+  void integrate (double t);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

--- a/liboctave/base-de.h
+++ b/liboctave/base-de.h
@@ -49,69 +49,11 @@
       return *this;
     }
 
-  // Derived classes must provide functions to actually do the
-  // integration.
-
-  // Return the vector of states at output time t.
-  virtual ColumnVector do_integrate (double tt) = 0;
+  // There must be a way for us to force the integration to restart.
 
-  // Return a matrix of states at each output time specified by t.
-  // The rows of the result matrix should each correspond to a new
-  // output time.
-  virtual Matrix do_integrate (const ColumnVector& tt) = 0;
-
-  virtual Matrix do_integrate (const ColumnVector& tt,
-			       const ColumnVector& ttcrit) = 0;
-
-  // There must also be a way for us to force the integration to
-  // restart.
   virtual void force_restart (void) = 0;
 
-  // Lots of ways to call the single function and optionally set and
-  // get additional information.
-
-  // Integrate to t from current point.
-  virtual ColumnVector integrate (double tt)
-    { return do_integrate (tt); }
-
-  // Set new x0, t0 and integrate to t.
-  virtual ColumnVector integrate (const ColumnVector& x0, double t0, double tt)
-    {
-      initialize (x0, t0);
-      return do_integrate (tt);
-    }
-
-  // Integrate from current point and return output at all points
-  // specified by t.
-  virtual Matrix integrate (const ColumnVector& tt)
-    { return do_integrate (tt); }
-
-  // Set new x0, t0 and integrate to return output at all points
-  // specified by t.
-  virtual Matrix integrate (const ColumnVector& x0, double t0,
-			    const ColumnVector& tt)
-    {
-      initialize (x0, t0);
-      return do_integrate (tt);
-    }
-
-  // Integrate from current point and return output at all points
-  // specified by t.
-  virtual Matrix integrate (const ColumnVector& tt,
-			    const ColumnVector& ttcrit)
-    { return do_integrate (tt, ttcrit); }
-
-  // Set new x0, t0 and integrate to return output at all points
-  // specified by t.
-  virtual Matrix integrate (const ColumnVector& x0, double t0,
-			    const ColumnVector& tt,
-			    const ColumnVector& ttcrit)
-    {
-      initialize (x0, t0);
-      return do_integrate (tt, ttcrit);
-    }
-
-  virtual void initialize (const ColumnVector& x0, double t0)
+  void initialize (const ColumnVector& x0, double t0)
     {
       x = x0;
       t = t0;

--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@
+2002-07-10  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* DLD-FUNCTIONS/odessa.cc: New file.
+	* Makefile.in (DLD_XSRC): Add it to the list.
+
 2002-07-05  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* pt-assign.cc (tree_multi_assignment::rvalue): Call

new file mode 100644
--- /dev/null
+++ b/src/DLD-FUNCTIONS/odessa.cc
@@ -0,0 +1,803 @@
+/*
+
+Copyright (C) 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave 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.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <string>
+
+#include <iomanip>
+#include <iostream>
+
+#include "ODESSA.h"
+#include "lo-mappers.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "ov-fcn.h"
+#include "pager.h"
+#include "pr-output.h"
+#include "unwind-prot.h"
+#include "utils.h"
+#include "variables.h"
+#include "parse.h"
+
+// Global pointer for user defined function required by odessa.
+static octave_function *odessa_f;
+static octave_function *odessa_j;
+static octave_function *odessa_b;
+
+static ODESSA_options odessa_opts;
+
+// Is this a recursive call?
+static int call_depth = 0;
+
+static ColumnVector
+odessa_user_f (double t, const ColumnVector& x, const ColumnVector& theta)
+{
+  ColumnVector retval;
+
+  octave_value_list args;
+
+  int n = x.length ();
+  int npar = theta.length ();
+
+  if (npar > 1)
+    args(2) = theta;
+  else if (npar == 1)
+    args(2) = theta(0);
+  else
+    args(2) = Matrix ();
+
+  if (n > 1)
+    args(1) = x;
+  else if (n == 1)
+    args(1) = x(0);
+  else
+    args(1) = Matrix ();
+
+  args(0) = t;
+
+  if (odessa_f)
+    {
+      octave_value_list tmp = odessa_f->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("odessa");
+	  return retval;
+	}
+
+      if (tmp.length () > 0 && tmp(0).is_defined ())
+	{
+	  retval = ColumnVector (tmp(0).vector_value ());
+
+	  if (error_state || retval.length () == 0)
+	    gripe_user_supplied_eval ("odessa");
+	}
+      else
+	gripe_user_supplied_eval ("odessa");
+    }
+
+  return retval;
+}
+
+static Matrix
+odessa_user_mf (double t, const ColumnVector& x, const ColumnVector& theta,
+		octave_function *mf)
+{
+  Matrix retval;
+
+  if (mf)
+    {
+      octave_value_list args;
+
+      int n = x.length ();
+      int npar = theta.length ();
+
+      if (npar > 1)
+	args(2) = theta;
+      else if (npar == 1)
+	args(2) = theta(0);
+      else
+	args(2) = Matrix ();
+
+      if (n > 1)
+	args(1) = x;
+      else if (n == 1)
+	args(1) = x(0);
+      else
+	args(1) = Matrix ();
+
+      args(0) = t;
+
+      octave_value_list tmp = mf->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("odessa");
+	  return retval;
+	}
+
+      if (tmp.length () > 0 && tmp(0).is_defined ())
+	{
+	  retval = tmp(0).matrix_value ();
+
+	  if (error_state || retval.length () == 0)
+	    gripe_user_supplied_eval ("odessa");
+	}
+      else
+	gripe_user_supplied_eval ("odessa");
+    }
+
+  return retval;
+}
+
+static Matrix
+odessa_user_j (double t, const ColumnVector& x, const ColumnVector& theta)
+{
+  return odessa_user_mf (t, x, theta, odessa_j);
+}
+
+static ColumnVector
+odessa_user_b (double t, const ColumnVector& x,
+	       const ColumnVector& theta, int column)
+{
+  ColumnVector retval;
+
+  if (odessa_b)
+    {
+      octave_value_list args;
+
+      int n = x.length ();
+      int npar = theta.length ();
+
+      args(3) = static_cast<double> (column);
+
+      if (npar > 1)
+	args(2) = theta;
+      else if (npar == 1)
+	args(2) = theta(0);
+      else
+	args(2) = Matrix ();
+
+      if (n > 1)
+	args(1) = x;
+      else if (n == 1)
+	args(1) = x(0);
+      else
+	args(1) = Matrix ();
+
+      args(0) = t;
+
+      octave_value_list tmp = odessa_b->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("odessa");
+	  return retval;
+	}
+
+      if (tmp.length () > 0 && tmp(0).is_defined ())
+	{
+	  retval = ColumnVector (tmp(0).vector_value ());
+
+	  if (error_state || retval.length () == 0)
+	    gripe_user_supplied_eval ("odessa");
+	}
+      else
+	gripe_user_supplied_eval ("odessa");
+    }
+
+  return retval;
+}
+
+static octave_value_list
+make_list (const Array<Matrix>& m_array)
+{
+  octave_value_list retval;
+
+  int len = m_array.length ();
+
+  retval.resize (len);
+
+  for (int i = 0; i < len; i++)
+    retval(i) = m_array(i);
+
+  return retval;
+}
+
+#define ODESSA_ABORT() \
+  do \
+    { \
+      unwind_protect::run_frame ("Fodessa"); \
+      return retval; \
+    } \
+  while (0)
+ 
+#define ODESSA_ABORT1(msg) \
+  do \
+    { \
+      ::error ("odessa: " msg); \
+      ODESSA_ABORT (); \
+    } \
+  while (0)
+
+#define ODESSA_ABORT2(fmt, arg) \
+  do \
+    { \
+      ::error ("odessa: " fmt, arg); \
+      ODESSA_ABORT (); \
+    } \
+  while (0)
+
+// --------------------------------
+// Everthing is so great above here
+// --------------------------------
+
+DEFUN_DLD (odessa, args, ,
+  "odessa (\"f\", x_0, theta, sx_0, t_out, t_crit)\n\
+\n\
+The string \"f\" may be substituted for the vector of strings\n\
+\n\
+               [\"f\"; \"j\"; \"b\"] \n\
+\n\
+You can use the function @code{odessa_options} to set optional\n\
+parameters for @code{odessa}.")
+{
+  octave_value_list retval;
+
+  unwind_protect::begin_frame ("Fodessa");
+
+  unwind_protect_int (call_depth);
+  call_depth++;
+
+  if (call_depth > 1)
+    ODESSA_ABORT1 ("invalid recursive call");
+
+  int nargin = args.length ();
+
+  if (nargin < 5 || nargin > 6)
+    {
+      print_usage ("odessa");
+      unwind_protect::run_frame ("Fodessa");
+      return retval;
+    }
+
+  odessa_f = 0;
+  odessa_j = 0;
+  odessa_b = 0;
+
+  octave_value f_arg = args(0);
+
+  if (f_arg.is_string ())
+    {
+      string_vector f_str_arg = f_arg.all_strings ();
+      
+      int len = f_str_arg.length ();
+      
+      if (len > 2)
+	{
+	  string t = f_str_arg(2);
+
+	  if (t.length () > 0)
+	    {
+	      odessa_b = is_valid_function (t, "odessa", true);
+	      
+	      if (! odessa_b)
+		ODESSA_ABORT1
+		  ("expecting b function name as argument 1");
+	    }
+	}
+
+      if (len > 1)
+	{
+	  string t = f_str_arg(1);
+	  
+	  if (t.length () > 0)
+	    {
+	      odessa_j = is_valid_function (t, "odessa", true);
+	      
+	      if (! odessa_j)
+		ODESSA_ABORT1
+		  ("expecting function name as argument 1");
+	    }
+	}
+      
+      if (len > 0)
+	odessa_f = is_valid_function (f_str_arg(0), "odessa", true);
+      
+      if (! odessa_f)
+	ODESSA_ABORT1 ("expecting function name as argument 1");
+    }
+  
+  ColumnVector state (args(1).vector_value ());
+
+  if (error_state)
+    ODESSA_ABORT1 ("expecting state vector as argument 2");
+
+  bool have_parameters = false;
+  
+  ColumnVector theta;
+  Matrix sensitivity_guess;
+
+  if (nargin == 5 || nargin == 6)
+    {
+      octave_value theta_arg = args(2);
+
+      if (! theta_arg.is_empty ())
+	{
+	  theta = ColumnVector (theta_arg.vector_value ());
+	  
+	  if (error_state)
+	    ODESSA_ABORT1
+	      ("expecting parameter vector as argument 3");
+	}
+
+      have_parameters = (theta.length () > 0);
+      
+      if (have_parameters)
+	{
+	  sensitivity_guess = args(3).matrix_value ();
+	  
+	  if (error_state)
+	    ODESSA_ABORT1
+	      ("expecting sensitivity guess as argument 4");
+	  
+	  if (sensitivity_guess.rows () != state.length ()
+	      || sensitivity_guess.columns () != theta.length ())
+	    ODESSA_ABORT1
+	      ("incorrect dimension for sensitivity guess");
+	}
+    }
+
+      ColumnVector out_times (args(4).vector_value ());
+
+      if (error_state)
+	ODESSA_ABORT1
+	  ("expecting output time vector as %s argument 5");
+
+      ColumnVector crit_times;
+
+      bool crit_times_set = false;
+
+      if (nargin==6)
+	{
+	  crit_times = ColumnVector (args(5).vector_value ());
+
+	  if (error_state)
+	    ODESSA_ABORT1
+	      ("expecting critical time vector as argument 6");
+
+	  crit_times_set = true;
+	}
+
+      ODESFunc func (odessa_user_f);
+
+      if (odessa_j)
+	func.set_jsub_function (odessa_user_j);
+
+      if (odessa_b)
+	func.set_bsub_function (odessa_user_b);
+
+      double tzero = out_times (0);
+
+      ODESSA_result output;
+
+      if (have_parameters)
+	{
+	  ODESSA dae = ODESSA (state, theta, sensitivity_guess, tzero, func);
+	  
+	  dae.copy (odessa_opts);
+
+	  if (crit_times_set)
+	    output = dae.integrate (out_times, crit_times);
+	  else
+	    output = dae.integrate (out_times);
+	}
+      else
+	{
+	  ODESSA dae = ODESSA (state, tzero, func);
+	  
+	  dae.copy (odessa_opts);
+	  
+	  if (crit_times_set)
+	    output = dae.integrate (out_times, crit_times);
+	  else
+	    output = dae.integrate (out_times);
+	}
+
+      if (! error_state)
+	{
+	  if (have_parameters)
+	    {
+	      retval(1) = make_list (output.state_sensitivity ());
+	    }
+	  
+	  retval(0) = output.state ();
+	}
+
+  unwind_protect::run_frame ("Fodessa");
+
+  return retval;
+}
+// -----------------------------
+// EVERYTHING SWELL BELOW HERE
+// -----------------------------
+
+
+typedef void (ODESSA_options::*da_set_opt_mf) (const Array<double>&);
+typedef void (ODESSA_options::*d_set_opt_mf) (double);
+typedef void (ODESSA_options::*i_set_opt_mf) (int);
+typedef void (ODESSA_options::*s_set_opt_mf) (const std::string&);
+
+typedef Array<double> (ODESSA_options::*da_get_opt_mf) (void) const;
+typedef double (ODESSA_options::*d_get_opt_mf) (void) const;
+typedef int (ODESSA_options::*i_get_opt_mf) (void) const;
+typedef std::string (ODESSA_options::*s_get_opt_mf) (void) const;
+
+#define MAX_TOKENS 3
+
+struct ODESSA_OPTIONS
+{
+  const char *keyword;
+  const char *kw_tok[MAX_TOKENS + 1];
+  int min_len[MAX_TOKENS + 1];
+  int min_toks_to_match;
+  da_set_opt_mf da_set_fcn;
+  d_set_opt_mf d_set_fcn;
+  i_set_opt_mf i_set_fcn;
+  s_set_opt_mf s_set_fcn;
+  da_get_opt_mf da_get_fcn;
+  d_get_opt_mf d_get_fcn;
+  i_get_opt_mf i_get_fcn;
+  s_get_opt_mf s_get_fcn;
+};
+
+static ODESSA_OPTIONS odessa_option_table [] =
+{
+  { "absolute tolerance",
+    { "absolute", "tolerance", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    &ODESSA_options::set_absolute_tolerance, 0, 0, 0,
+    &ODESSA_options::absolute_tolerance, 0, 0, 0, },
+
+  { "initial step size",
+    { "initial", "step", "size", 0, },
+    { 3, 0, 0, 0, }, 1,
+    0, &ODESSA_options::set_initial_step_size, 0, 0,
+    0, &ODESSA_options::initial_step_size, 0, 0, },
+
+  { "integration method",
+    { "integration", "method", 0, 0, },
+    { 3, 0, 0, 0, }, 1,
+    0, 0, 0, &ODESSA_options::set_integration_method,
+    0, 0, 0, &ODESSA_options::integration_method, },
+
+  { "maximum step size",
+    { "maximum", "step", "size", 0, },
+    { 2, 0, 0, 0, }, 1,
+    0, &ODESSA_options::set_maximum_step_size, 0, 0,
+    0, &ODESSA_options::maximum_step_size, 0, 0, },
+
+  { "minimum step size",
+    { "minimum", "step", "size", 0, },
+    { 2, 0, 0, 0, }, 1,
+    0, &ODESSA_options::set_minimum_step_size, 0, 0,
+    0, &ODESSA_options::minimum_step_size, 0, 0, },
+
+  { "relative tolerance",
+    { "relative", "tolerance", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    0, &ODESSA_options::set_relative_tolerance, 0, 0,
+    0, &ODESSA_options::relative_tolerance, 0, 0, },
+
+  { "step limit",
+    { "step", "limit", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    0, 0, &ODESSA_options::set_step_limit, 0,
+    0, 0, &ODESSA_options::step_limit, 0, },
+
+  { 0,
+    { 0, 0, 0, 0, },
+    { 0, 0, 0, 0, }, 0,
+    0, 0, 0, 0,
+    0, 0, 0, 0, },
+};
+
+static void
+print_odessa_option_list (std::ostream& os)
+{
+  print_usage ("odessa_options", 1);
+
+  os << "\n"
+     << "Options for odessa include:\n\n"
+     << "  keyword                                  value\n"
+     << "  -------                                  -----\n\n";
+
+  ODESSA_OPTIONS *list = odessa_option_table;
+
+  const char *keyword;
+  while ((keyword = list->keyword) != 0)
+    {
+      os << "  "
+	 << std::setiosflags (std::ios::left) << std::setw (40)
+	 << keyword
+	 << std::resetiosflags (std::ios::left)
+	 << " ";
+
+      if (list->da_get_fcn)
+	{
+	  Array<double> val = (odessa_opts.*list->da_get_fcn) ();
+	  if (val.length () == 1)
+	    {
+	      if (val(0) < 0.0)
+		os << "computed automatically";
+	      else
+		os << val(0);
+	    }
+	  else
+	    {
+	      os << "\n\n";
+	      // XXX FIXME XXX
+	      Matrix tmp = Matrix (ColumnVector (val));
+	      octave_print_internal (os, tmp, false, 2);
+	      os << "\n";
+	    }
+	}
+      else if (list->d_get_fcn)
+	{
+	  double val = (odessa_opts.*list->d_get_fcn) ();
+	  if (val < 0.0)
+	    os << "computed automatically";
+	  else
+	    os << val;
+	}
+      else if (list->i_get_fcn)
+	{
+	  int val = (odessa_opts.*list->i_get_fcn) ();
+	  if (val < 0)
+	    os << "infinite";
+	  else
+	    os << val;
+	}
+      else if (list->s_get_fcn)
+	{
+	  os << (odessa_opts.*list->s_get_fcn) ();
+	}
+      else
+	panic_impossible ();
+
+      os << "\n";
+      list++;
+    }
+
+  os << "\n";
+}
+
+static void
+set_odessa_option (const std::string& keyword, double val)
+{
+  ODESSA_OPTIONS *list = odessa_option_table;
+
+  while (list->keyword != 0)
+    {
+      if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
+				list->min_toks_to_match, MAX_TOKENS))
+	{
+	  if (list->da_set_fcn)
+	    {
+	      Array<double> tmp (1, val);
+	      (odessa_opts.*list->da_set_fcn) (tmp);
+	    }
+	  else if (list->d_set_fcn)
+	    {
+	      (odessa_opts.*list->d_set_fcn) (val);
+	    }
+	  else
+	    {
+	      if (xisnan (val))
+		{
+		  error ("odessa_options: %s: expecting integer, found NaN",
+			 keyword.c_str ());
+		}
+	      else
+		(odessa_opts.*list->i_set_fcn) (NINT (val));
+	    }
+	  return;
+	}
+      list++;
+    }
+
+  warning ("odessa_options: no match for `%s'", keyword.c_str ());
+}
+
+static void
+set_odessa_option (const std::string& keyword, const Array<double>& val)
+{
+  ODESSA_OPTIONS *list = odessa_option_table;
+
+  while (list->keyword != 0)
+    {
+      if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
+				list->min_toks_to_match, MAX_TOKENS))
+	{
+	  if (list->da_set_fcn)
+	    (odessa_opts.*list->da_set_fcn) (val);
+	  else
+	    error ("odessa_options: no function to handle vector option");
+
+	  return;
+	}
+      list++;
+    }
+
+  warning ("odessa_options: no match for `%s'", keyword.c_str ());
+}
+
+static void
+set_odessa_option (const std::string& keyword, const std::string& val)
+{
+  ODESSA_OPTIONS *list = odessa_option_table;
+
+  while (list->keyword != 0)
+    {
+      if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
+				list->min_toks_to_match, MAX_TOKENS))
+	{
+	  if (list->s_set_fcn)
+	    (odessa_opts.*list->s_set_fcn) (val);
+	  else
+	    error ("odessa_options: no function to handle string option");
+
+	  return;
+	}
+      list++;
+    }
+
+  warning ("odessa_options: no match for `%s'", keyword.c_str ());
+}
+
+static octave_value_list
+show_odessa_option (const std::string& keyword)
+{
+  octave_value retval;
+
+  ODESSA_OPTIONS *list = odessa_option_table;
+
+  while (list->keyword != 0)
+    {
+      if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
+				list->min_toks_to_match, MAX_TOKENS))
+	{
+	  if (list->da_get_fcn)
+	    {
+	      Array<double> val = (odessa_opts.*list->da_get_fcn) ();
+	      if (val.length () == 1)
+		{
+		  if (val(0) < 0.0)
+		    retval = "computed automatically";
+		  else
+		    retval = val(0);
+		}
+	      else
+		retval = ColumnVector (val);
+	    }
+	  else if (list->d_get_fcn)
+	    {
+	      double val = (odessa_opts.*list->d_get_fcn) ();
+	      if (val < 0.0)
+		retval = "computed automatically";
+	      else
+		retval = val;
+	    }
+	  else if (list->i_get_fcn)
+	    {
+	      int val = (odessa_opts.*list->i_get_fcn) ();
+	      if (val < 0)
+		retval = "infinite";
+	      else
+		retval = static_cast<double> (val);
+	    }
+	  else if (list->s_get_fcn)
+	    {
+	      retval = (odessa_opts.*list->s_get_fcn) ();
+	    }
+	  else
+	    panic_impossible ();
+
+	  return retval;
+	}
+      list++;
+    }
+
+  warning ("odessa_options: no match for `%s'", keyword.c_str ());
+
+  return retval;
+}
+
+DEFUN_DLD (odessa_options, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} odessa_options (@var{opt}, @var{val})\n\
+When called with two arguments, this function allows you set options\n\
+parameters for the function @code{odessa}.  Given one argument,\n\
+@code{odessa_options} returns the value of the corresponding option.  If\n\
+no arguments are supplied, the names of all the available options and\n\
+their current values are displayed.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+
+  if (nargin == 0)
+    {
+      print_odessa_option_list (octave_stdout);
+    }
+  else if (nargin == 1 || nargin == 2)
+    {
+      std::string keyword = args(0).string_value ();
+
+      if (! error_state)
+	{
+	  if (nargin == 1)
+	    retval = show_odessa_option (keyword);
+	  else
+	    {
+	      if (args(1).is_string ())
+		{
+		  std::string val = args(1).string_value ();
+
+		  if (! error_state)
+		    set_odessa_option (keyword, val);
+		}
+	      else if (args(1).is_scalar_type ())
+		{
+		  double val = args(1).double_value ();
+
+		  if (! error_state)
+		    set_odessa_option (keyword, val);
+		}
+	      else
+		{
+		  Array<double> val = args(1).vector_value ();
+
+		  if (! error_state)
+		    set_odessa_option (keyword, val);
+		}
+	    }
+	}
+    }
+  else
+    print_usage ("odessa_options");
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -39,12 +39,13 @@
   endif
 endif
 
-DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc daspk.cc \
-	dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc filter.cc find.cc \
-	fsolve.cc gammainc.cc getgrent.cc getpwent.cc getrusage.cc \
-	givens.cc hess.cc ifft.cc ifft2.cc inv.cc kron.cc log.cc lpsolve.cc \
-	lsode.cc lu.cc minmax.cc pinv.cc qr.cc quad.cc qz.cc rand.cc \
-	schur.cc sort.cc svd.cc syl.cc time.cc
+DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc \
+	daspk.cc dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc \
+	filter.cc find.cc fsolve.cc gammainc.cc getgrent.cc \
+	getpwent.cc getrusage.cc givens.cc hess.cc ifft.cc ifft2.cc \
+	inv.cc kron.cc log.cc lpsolve.cc lsode.cc lu.cc minmax.cc \
+	odessa.cc pinv.cc qr.cc quad.cc qz.cc rand.cc schur.cc \
+	sort.cc svd.cc syl.cc time.cc
 
 DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))