--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,15 @@
+2002-07-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* DAE.cc: Delete.
+
+	* DAERT.h, DAERTFunc.h, DASRT.h, DASRT.cc: New files for DAE
+	solving with root finding.
+	* Makefile.in: Add them to the appropriate lists.
+
+	* base-dae.h: New file.
+	* Makefile.in (INCLUDES): Add it to the list.
+	* DAE.h (DAE): Derive from base_diff_alg_eqn, not base_diff_eqn.
+
 2002-07-10  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* ODE.h: Move integrate and do_integrate method declarations and

deleted file mode 100644
--- a/liboctave/DAE.cc
+++ /dev/null
@@ -1,68 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 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 "DAE.h"
-#include "lo-error.h"
-
-DAE::DAE (const ColumnVector& xx, const ColumnVector& xxdot,
-	  double tt, DAEFunc& f)
-  : base_diff_eqn (xx, tt), DAEFunc (f), xdot (xxdot)
-{
-  if (x.length () != xdot.length ())
-    ; // XXX FIXME XXX -- exception!
-}
-
-void
-DAE::initialize (const ColumnVector& xx, double tt)
-{
-  if (xx.length () != xdot.length ())
-    ; // XXX FIXME XXX -- exception!
-  else
-    base_diff_eqn::initialize (xx, tt);
-}
-
-void
-DAE::initialize (const ColumnVector& xx, const ColumnVector& xxdot,
-		 double tt)
-{
-  if (xx.length () != xxdot.length ())
-    ; // XXX FIXME XXX -- exception!
-  else
-    {
-      base_diff_eqn::initialize (xx, tt);
-      xdot = xxdot;
-    }
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/

--- a/liboctave/DAE.h
+++ b/liboctave/DAE.h
@@ -23,54 +23,38 @@
 #if !defined (octave_DAE_h)
 #define octave_DAE_h 1
 
-#if defined (__GNUG__)
-#pragma interface
-#endif
-
 #include "DAEFunc.h"
-#include "base-de.h"
+#include "base-dae.h"
 
 class
-DAE : public base_diff_eqn, public DAEFunc
+DAE : public base_diff_alg_eqn, public DAEFunc
 {
 public:
 
   DAE (void)
-    : base_diff_eqn (), DAEFunc (), xdot () { }
+    : base_diff_alg_eqn (), DAEFunc () { }
 
   DAE (const ColumnVector& xx, double tt, DAEFunc& f)
-    : base_diff_eqn (xx, tt), DAEFunc (f), xdot (x.capacity (), 0.0) { }
+    : base_diff_alg_eqn (xx, tt), DAEFunc (f) { }
 
   DAE (const ColumnVector& xx, const ColumnVector& xxdot,
-       double tt, DAEFunc& f);
+       double tt, DAEFunc& f)
+    : base_diff_alg_eqn (xx, xxdot, tt), DAEFunc (f) { }
 
   DAE (const DAE& a)
-    : base_diff_eqn (a), DAEFunc (a), xdot (a.xdot) { }
+    : base_diff_alg_eqn (a), DAEFunc (a){ }
 
   DAE& operator = (const DAE& a)
     {
       if (this != &a)
 	{
-	  base_diff_eqn::operator = (a);
+	  base_diff_alg_eqn::operator = (a);
 	  DAEFunc::operator = (a);
-
-	  xdot = a.xdot;
 	}
       return *this;
     }
 
   ~DAE (void) { }
-
-  ColumnVector state_derivative (void) { return xdot; }
-
-  void initialize (const ColumnVector& xx, double tt);
-
-  void initialize (const ColumnVector& xx, const ColumnVector& xxdot,
-		   double tt);
-
-protected:
-
-  ColumnVector xdot;
 };
 
 #endif

--- a/liboctave/DAEFunc.h
+++ b/liboctave/DAEFunc.h
@@ -39,10 +39,10 @@
 
   typedef ColumnVector (*DAERHSFunc) (const ColumnVector& x,
 				      const ColumnVector& xdot,
-				      double, int&); 
+				      double t, int& ires);
 
   typedef DAEJac (*DAEJacFunc) (const ColumnVector& x,
-				const ColumnVector& xdot, double);
+				const ColumnVector& xdot, double t);
 
   DAEFunc (void)
     : fun (0), jac (0) { }

new file mode 100644
--- /dev/null
+++ b/liboctave/DAERT.h
@@ -0,0 +1,70 @@
+/*
+
+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_DAERT_h)
+#define octave_DAERT_h 1
+
+#include "DAE.h"
+#include "DAERTFunc.h"
+#include "base-dae.h"
+
+class
+DAERT : public base_diff_alg_eqn, public DAERTFunc
+{
+public:
+
+  DAERT (void)
+    : base_diff_alg_eqn (), DAERTFunc () { }
+
+  DAERT (const ColumnVector& x, const ColumnVector& xdot, double t,
+	DAERTFunc& f)
+    : base_diff_alg_eqn (x, xdot, t), DAERTFunc (f) { }
+
+  DAERT (const DAERT& a)
+    : base_diff_alg_eqn (a), DAERTFunc (a) { }
+
+  DAERT& operator = (const DAERT& a)
+    {
+      if (this != &a)
+	{
+	  base_diff_alg_eqn::operator = (a);
+	  DAERTFunc::operator = (a);
+
+	}
+      return *this;
+    }
+
+  ~DAERT (void) { }
+
+  void initialize (const ColumnVector& x, const ColumnVector& xdot, double t)
+    {
+      base_diff_alg_eqn::initialize (x, xdot, t);
+    }
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/DAERTFunc.h
@@ -0,0 +1,84 @@
+/*
+
+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_DAERTFunc_h)
+#define octave_DAERTFunc_h 1
+
+#include "dMatrix.h"
+
+class
+DAERTFunc : DAEFunc
+{
+public:
+
+  typedef ColumnVector (*DAERTConstrFunc) (const ColumnVector& x, double t);
+
+  DAERTFunc (void)
+    : DAEFunc (), constr (0) { }
+
+  DAERTFunc (DAERHSFunc f)
+    : DAEFunc (f), constr (0) { }
+
+  DAERTFunc (DAERHSFunc f, DAEJacFunc j)
+    : DAEFunc (f, j), constr (0) { }
+
+  DAERTFunc (DAERHSFunc f, DAERTConstrFunc cf)
+    : DAEFunc (f), constr (cf) { }
+
+  DAERTFunc (DAERHSFunc f, DAERTConstrFunc cf, DAEJacFunc j)
+    : DAEFunc (f, j), constr (cf) { }
+
+  DAERTFunc (const DAERTFunc& a)
+    : DAEFunc (a), constr (a.constr) { }
+
+  DAERTFunc& operator = (const DAERTFunc& a)
+    {
+      if (this != &a)
+	{
+	  DAEFunc::operator = (a);
+	  constr = a.constr;
+	}
+      return *this;
+    }
+
+  ~DAERTFunc (void) { }
+
+  DAERTConstrFunc constraint_function (void) const { return constr; }
+
+  DAERTFunc& set_constraint_function (DAERTConstrFunc cf)
+    {
+      constr = cf;
+      return *this;
+    }
+
+protected:
+
+  DAERTConstrFunc constr;
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

--- a/liboctave/DASPK.cc
+++ b/liboctave/DASPK.cc
@@ -195,7 +195,7 @@
 
   // Fix up the matrix of partial derivatives for daspk.
 
-  tmp_dfdx = tmp_dfdx +  cj * tmp_dfdxdot;
+  tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot;
 
   for (int j = 0; j < nn; j++)
     for (int i = 0; i < nn; i++)

new file mode 100644
--- /dev/null
+++ b/liboctave/DASRT.cc
@@ -0,0 +1,667 @@
+/*
+
+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 <fstream.h>
+
+#include <cstdlib>
+#include <cfloat>
+#include <cmath>
+#include "defun-dld.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"
+
+// For instantiating the Array<Matrix> object.
+#include "Array.h"
+#include "Array.cc"
+
+#include "DASRT.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+
+#ifndef F77_FUNC
+#define F77_FUNC(x, X) F77_FCN (x, X)
+#endif
+
+extern "C"
+{
+  int F77_FUNC (ddasrt, DASRT) (int (*)(const double&, double*, double*,
+					double*, int&, double*, int*),
+				const int&, const double&, double*, double*,
+				const double&, int*, double*, double*,
+				int&, double*, const int&, int*, 
+				const int&, double*, int*,
+				int (*)(const double&, double*,
+					double*, double*,
+					const double&, double*, int*),
+				int (*)(const int&, const double&, double*,
+					const int&, double*, double*, int*),
+				const int&, int*);
+}
+
+template class Array<Matrix>;
+
+static DAEFunc::DAERHSFunc user_fsub;
+static DAEFunc::DAEJacFunc user_jsub;
+static DAERTFunc::DAERTConstrFunc user_csub;
+static int nn;
+
+static int
+ddasrt_f (const double& t, double *state, double *deriv, double *delta,
+          int& ires, double *rpar, int *ipar)
+{
+  ColumnVector tmp_state (nn);
+  for (int i = 0; i < nn; i++)
+    tmp_state(i) = state[i];
+
+  ColumnVector tmp_deriv (nn);
+  for (int i = 0; i < nn; i++)
+    tmp_deriv(i) = deriv[i];
+
+  ColumnVector tmp_fval = user_fsub (tmp_state, tmp_deriv, t, ires);
+
+  if (tmp_fval.length () == 0)
+    ires = -2;
+  else
+    {
+      for (int i = 0; i < nn; i++)
+	delta[i] = tmp_fval(i);
+    }
+
+  return 0;
+}
+
+//typedef int (*efptr) (const double& t, const int& n, double *state,
+//		      double *ework, double *rpar, int *ipar,
+//		      const int& ieform, int& ires);
+
+//static efptr e_fun;
+
+static int
+ddasrt_j (const double& t, double *state, double *deriv,
+	  double *pdwork, const double& cj, double *rpar, int *ipar) 
+{
+  ColumnVector tmp_state (nn);
+  for (int i = 0; i < nn; i++)
+    tmp_state(i) = state[i];
+
+  ColumnVector tmp_deriv (nn);
+  for (int i = 0; i < nn; i++)
+    tmp_deriv(i) = deriv[i];
+
+  // XXX FIXME XXX
+
+  Matrix tmp_dfdxdot (nn, nn);
+  Matrix tmp_dfdx (nn, nn);
+
+  DAEFunc::DAEJac tmp_jac;
+  tmp_jac.dfdxdot = &tmp_dfdxdot;
+  tmp_jac.dfdx    = &tmp_dfdx;
+
+  tmp_jac = user_jsub (tmp_state, tmp_deriv, t);
+
+  // Fix up the matrix of partial derivatives for dasrt.
+
+  tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot;
+
+  for (int j = 0; j < nn; j++)
+    for (int i = 0; i < nn; i++)
+      pdwork[j*nn+i] = tmp_dfdx.elem (i, j);
+
+  return 0;
+}
+
+static int
+ddasrt_g (const int& neq, const double& t, double *state, const int& ng,
+	  double *gout, double *rpar, int *ipar) 
+{
+  int n = neq;
+
+  ColumnVector tmp_state (n);
+  for (int i = 0; i < n; i++)
+    tmp_state(i) = state[i];
+
+  ColumnVector tmp_fval = user_csub (tmp_state, t);
+
+  for (int i = 0; i < ng; i++)
+    gout[i] = tmp_fval(i);
+
+  return 0;
+}
+
+
+DASRT::DASRT (void)
+  : DAERT ()
+{
+  initialized = false;
+  restart = false;
+
+  stop_time_set = false;
+  stop_time = 0.0;
+
+  sanity_checked = false;
+
+  info.resize (30, 0);
+
+  npar = 0;
+
+  liw = 0;
+  lrw = 0;
+}
+
+DASRT::DASRT (const int& ng, const ColumnVector& state, 
+	      const ColumnVector& deriv, double time, DAERTFunc& f)
+  : DAERT (state, deriv, time, f)
+{
+  n = size ();
+
+  initialized = false;
+  restart = false;
+
+  stop_time_set = false;
+  stop_time = 0.0;
+
+  DAERTFunc::operator = (f);
+
+  sanity_checked = false;
+
+  info.resize (30, 0);
+  jroot.resize (ng, 1);
+
+  npar = 0;
+
+  rpar.resize (npar+1);
+  ipar.resize (npar+1);
+
+  info(11) = npar;
+
+  // Also store it here, for communication with user-supplied
+  // subroutines.
+  ipar(0) = npar;
+
+  y.resize (n, 1, 0.0);
+  ydot.resize (n, 1, 0.0);
+}
+
+void
+DASRT::init_work_size (int info_zero)
+{
+  double t;
+  double *py = y.fortran_vec ();
+  double *pydot = ydot.fortran_vec ();
+  double rel_tol = relative_tolerance ();
+  double abs_tol = absolute_tolerance ();
+  int *pinfo = info.fortran_vec ();
+  double *prpar = rpar.fortran_vec ();
+  int *pipar = ipar.fortran_vec ();
+  int *pjroot = jroot.fortran_vec ();
+  int idid;
+
+  // We do not have to lie.
+  rwork.resize (5000+9*n+n*n, 0.0);
+  iwork.resize (n+20, 0);
+
+  liw = n+20;
+  lrw = 5000+9*n+n*n;
+
+  double *prwork = rwork.fortran_vec ();
+  int *piwork = iwork.fortran_vec ();
+
+
+  F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo,
+			    &rel_tol, &abs_tol, idid, prwork, lrw,
+			    piwork, liw, prpar, pipar, ddasrt_j,
+			    ddasrt_g, ng, pjroot);
+
+  int iwadd = iwork(18);
+
+
+  if (iwadd > 0)
+    liw += iwadd;
+
+  info(0) = 0;
+
+  iwork.resize (liw, 0);
+
+  piwork = iwork.fortran_vec ();
+
+  F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo,
+			    &rel_tol, &abs_tol, idid, prwork, lrw,
+			    piwork, liw, prpar, pipar, ddasrt_j,
+			    ddasrt_g, ng, pjroot);
+
+  int rwadd = iwork(19);
+
+
+  if (rwadd > 0)
+    lrw += rwadd;
+
+  rwork.resize (lrw, 0.0);
+
+  info(0) = info_zero;
+
+}
+
+void
+DASRT::force_restart (void)
+{
+  restart = true;
+  integration_error = false;
+}
+
+void
+DASRT::set_stop_time (double t)
+{
+  stop_time_set = true;
+  stop_time = t;
+}
+
+void
+DASRT::set_ng (int the_ng)
+{
+  ng = the_ng;
+}
+
+int DASRT::get_ng (void)
+{
+  return ng;
+}
+
+void
+DASRT::clear_stop_time (void)
+{
+  stop_time_set = false;
+}
+
+void
+DASRT::integrate (double tout)
+{
+  DASRT_result retval;
+
+  if (! initialized)
+    {
+      info(0) = 0;
+
+      for (int i = 0; i < n; i++)
+       {
+	y(i,0) = x(i);
+        ydot(i,0) = xdot(i);
+       }
+
+      integration_error = false;
+
+      user_fsub = DAEFunc::function ();
+      user_jsub = DAEFunc::jacobian_function ();
+      user_csub = DAERTFunc::constraint_function ();
+
+      if (user_jsub)
+	info(4) = 1;
+      else
+	info(4) = 0;
+
+      if (! sanity_checked)
+	{
+	  int ires = 0;
+
+	  ColumnVector fval = user_fsub (x, xdot, t, ires);
+
+	  if (fval.length () != x.length ())
+	    {
+	      (*current_liboctave_error_handler)
+		("dassl: inconsistent sizes for state and residual vectors");
+
+	      integration_error = true;
+	      return;
+	    }
+
+	  sanity_checked = true;
+	}
+
+  
+      init_work_size (info(0));
+
+
+
+
+      if (iwork.length () != liw)
+	iwork.resize (liw);
+
+      if (rwork.length () != lrw)
+	rwork.resize (lrw);
+
+      abs_tol = absolute_tolerance ();
+      rel_tol = relative_tolerance ();
+
+
+      if (initial_step_size () >= 0.0)
+	{
+	  rwork(2) = initial_step_size ();
+	  info(7) = 1;
+	}
+      else
+	info(7) = 0;
+
+      if (step_limit () >= 0)
+	{
+	  info(11) = 1;
+	  iwork(18) = step_limit ();
+	}
+      else
+	info(11) = 0;
+
+      if (maximum_step_size () >= 0.0)
+	{
+	  rwork(1) = maximum_step_size ();
+	  info(6) = 1;
+	}
+      else
+	info(6) = 0;
+
+
+      py = y.fortran_vec ();
+      pydot = ydot.fortran_vec ();
+      pinfo = info.fortran_vec ();
+      piwork = iwork.fortran_vec ();
+      prwork = rwork.fortran_vec ();
+      prpar = rpar.fortran_vec ();
+      pipar = ipar.fortran_vec ();
+      pjroot = jroot.fortran_vec ();
+
+      info(5) = 0;
+      info(8) = 0;
+      initialized = true;
+    }
+
+  if (restart)
+    {
+      info(0) = 0;
+
+      if (stop_time_set)
+	{
+	  info(3) = 1;
+	  rwork(0) = stop_time;
+	}
+      else
+	info(3) = 0;
+    }
+
+
+
+
+  F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, py, pydot, tout, pinfo,
+			    &rel_tol, &abs_tol, idid, prwork, lrw,
+			    piwork, liw, prpar, pipar, ddasrt_j,
+			    ddasrt_g, ng, pjroot));
+
+  if (f77_exception_encountered)
+    {
+      integration_error = true;
+      (*current_liboctave_error_handler) ("unrecoverable error in dassl");
+    }
+  else
+    {
+      switch (idid)
+	{
+	case 0: // Initial conditions made consistent.
+	case 1: // A step was successfully taken in intermediate-output
+	        // mode. The code has not yet reached TOUT.
+	case 2: // The integration to TOUT was successfully completed
+	        // (T=TOUT) by stepping exactly to TOUT.
+	case 3: // The integration to TOUT was successfully completed
+	        // (T=TOUT) by stepping past TOUT.  Y(*) is obtained by
+	        // interpolation.  YPRIME(*) is obtained by interpolation.
+	case 5: // The integration to TSTOP was successfully completed
+	        // (T=TSTOP) by stepping to TSTOP within the
+	        // tolerance.  Must restart to continue.
+	  for (int i = 0; i < n; i++)
+	    x(i) = y(i,0);
+	  t = tout;
+	  break;
+
+	case 4: //  We've hit the stopping condition.
+          for (int i = 0; i < n; i++)
+            x(i) = y(i,0);
+          break;
+
+	case -1: // A large amount of work has been expended.  (~500 steps).
+	case -2: // The error tolerances are too stringent.
+	case -3: // The local error test cannot be satisfied because you
+	         // specified a zero component in ATOL and the
+		 // corresponding computed solution component is zero.
+		 // Thus, a pure relative error test is impossible for
+		 // this component.
+	case -6: // DDASRT had repeated error test failures on the last
+		 // attempted step.
+	case -7: // The corrector could not converge.
+	case -8: // The matrix of partial derivatives is singular.
+	case -9: // The corrector could not converge.  There were repeated
+		 // error test failures in this step.
+	case -10: // The corrector could not converge because IRES was
+		  // equal to minus one.
+	case -11: // IRES equal to -2 was encountered and control is being
+		  // returned to the calling program.
+	case -12: // DASSL failed to compute the initial YPRIME.
+	case -33: // The code has encountered trouble from which it cannot
+		  // recover. A message is printed explaining the trouble
+		  // and control is returned to the calling program. For
+		  // example, this occurs when invalid input is detected.
+	default:
+	  integration_error = true;
+	  (*current_liboctave_error_handler)
+	    ("ddasrt failed with IDID = %d", idid);
+	  break;
+	}
+    }
+}
+
+DASRT_result
+DASRT::integrate (const ColumnVector& tout)
+{
+  DASRT_result retval;
+
+  Matrix x_out;
+  Matrix xdot_out;
+  ColumnVector t_out;
+
+  int oldj = 0;
+
+  int n_out = tout.capacity ();
+
+
+  if (n_out > 0 && n > 0)
+    {
+      x_out.resize (n_out, n);
+      xdot_out.resize (n_out, n);
+      t_out.resize (n_out);
+
+      for (int j = 0; j < n_out; j++)
+	{
+	  integrate (tout(j));
+	  if (integration_error)
+	    {
+	      retval = DASRT_result (x_out, xdot_out, t_out);
+	      return retval;
+	    }
+          if (idid == 4)
+            t_out(j) = t;
+          else
+            t_out(j) = tout(j);
+
+
+	  for (int i = 0; i < n; i++)
+	    {
+	      x_out(j,i) = y(i,0);
+	      xdot_out(j,i) = ydot(i,0);
+	    }
+          if (idid ==4)
+           {
+            oldj = j;
+            j = n_out;
+            x_out.resize (oldj+1, n);
+            xdot_out.resize (oldj+1, n);
+            t_out.resize (oldj+1);
+           }
+	}
+    }
+
+  retval = DASRT_result (x_out, xdot_out, t_out);
+
+  return retval;
+}
+
+DASRT_result
+DASRT::integrate (const ColumnVector& tout, const ColumnVector& tcrit) 
+{
+  DASRT_result retval;
+
+  Matrix x_out;
+  Matrix xdot_out;
+  ColumnVector t_outs;
+
+  int n_out = tout.capacity ();
+
+  if (n_out > 0 && n > 0)
+    {
+      x_out.resize (n_out, n);
+      xdot_out.resize (n_out, n);
+      t_outs.resize (n_out);
+
+      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 = 1;
+		  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 = 0;
+		      i_crit++;
+		      do_restart = true;
+		    }
+		  else
+		    {
+		      clear_stop_time ();
+		      t_out = next_out;
+		      save_output = 1;
+		      i_out++;
+		    }
+		}
+	      else
+		{
+		  set_stop_time (next_crit);
+		  t_out = next_out;
+		  save_output = 1;
+		  i_out++;
+		}
+
+	      integrate (t_out);
+
+	      if (integration_error)
+		{
+		  retval = DASRT_result (x_out, xdot_out, t_outs);
+		  return retval;
+		}
+              if (idid == 4)
+                t_out = t;
+
+	      if (save_output)
+		{
+		  for (int i = 0; i < n; i++)
+		    {
+		      x_out(i_out-1,i) = y(i,0);
+		      xdot_out(i_out-1,i) = ydot(i,0);
+		    }
+                  t_outs(i_out-1) = t_out;
+                  if (idid ==4)
+                    {
+                      x_out.resize (i_out, n);
+                      xdot_out.resize (i_out, n);
+                      t_outs.resize (i_out);
+                      i_out = n_out;
+                    }
+
+		}
+
+	      if (do_restart)
+		force_restart ();
+	    }
+
+	  retval = DASRT_result (x_out, xdot_out, t_outs);
+	}
+      else
+	{
+	  retval = integrate (tout);
+
+	  if (integration_error)
+	    return retval;
+	}
+    }
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

new file mode 100644
--- /dev/null
+++ b/liboctave/DASRT.h
@@ -0,0 +1,235 @@
+/*
+
+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_DASRT_h)
+#define octave_DASRT_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#include <cfloat>
+#include <cmath>
+
+#include "DAERT.h"
+
+class
+DASRT_options
+{
+public:
+
+  DASRT_options (void) { init (); }
+
+  DASRT_options (const DASRT_options& opt) { copy (opt); }
+
+  DASRT_options& operator = (const DASRT_options& opt)
+    {
+      if (this != &opt)
+	copy (opt);
+
+      return *this;
+    }
+
+  ~DASRT_options (void) { }
+
+  void init (void)
+    {
+      double sqrt_eps = ::sqrt (DBL_EPSILON);
+      x_absolute_tolerance = sqrt_eps;
+      x_initial_step_size = -1.0;
+      x_maximum_step_size = -1.0;
+      x_minimum_step_size = 0.0;
+      x_relative_tolerance = sqrt_eps;
+      x_step_limit = -1;
+    }
+
+  void copy (const DASRT_options& opt)
+    {
+      x_absolute_tolerance = opt.x_absolute_tolerance;
+      x_initial_step_size = opt.x_initial_step_size;
+      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 = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); }
+
+  void set_initial_step_size (double val)
+    { x_initial_step_size = (val >= 0.0) ? val : -1.0; }
+
+  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 >= 0) ? val : -1; }
+
+  double absolute_tolerance (void) { return x_absolute_tolerance; }
+
+  double initial_step_size (void) { return x_initial_step_size; }
+
+  double maximum_step_size (void) { return x_maximum_step_size; }
+
+  double minimum_step_size (void) { return x_minimum_step_size; }
+
+  double relative_tolerance (void) { return x_relative_tolerance; }
+
+  int step_limit (void) { return x_step_limit; }
+
+private:
+
+  double x_absolute_tolerance;
+  double x_initial_step_size;
+  double x_maximum_step_size;
+  double x_minimum_step_size;
+  double x_relative_tolerance;
+  int x_step_limit;
+};
+
+class
+DASRT_result
+{
+public:
+
+  DASRT_result (void) { }
+
+  DASRT_result (const Matrix& xx, const Matrix& xxdot, const ColumnVector& tt)
+    : x (xx), xdot (xxdot), t (tt) { }
+
+  DASRT_result (const DASRT_result& r)
+    : x (r.x), xdot (r.xdot), t (r.t) { }
+
+  DASRT_result& operator = (const DASRT_result& r)
+    {
+      if (this != &r)
+	{
+	  x = r.x;
+	  xdot = r.xdot;
+          t = r.t;
+	}
+      return *this;
+    }
+
+  ~DASRT_result (void) { }
+
+  Matrix state (void) const { return x; }
+  Matrix deriv (void) const { return xdot; }
+  ColumnVector times (void) const { return t; }
+
+private:
+
+  Matrix x;
+  Matrix xdot;
+  ColumnVector t;
+};
+
+class
+DASRT : public DAERT, public DASRT_options
+{
+public:
+
+  DASRT (void);
+
+  DASRT (const int& ng, const ColumnVector& x, const ColumnVector& xdot, 
+	 double time, DAERTFunc& f);
+
+  ~DASRT (void) { }
+
+  void force_restart (void);
+
+  void set_stop_time (double t);
+  void clear_stop_time (void);
+  void set_ng (int the_ng);
+  int get_ng (void);
+
+  DASRT_result integrate (const ColumnVector& tout);
+
+  DASRT_result integrate (const ColumnVector& tout,
+			  const ColumnVector& tcrit); 
+
+private:
+
+  bool initialized;
+
+  bool sanity_checked;
+
+  bool stop_time_set;
+  double stop_time;
+
+  bool restart;
+
+  bool integration_error;
+
+  int liw;  
+  int lrw;
+  int idid;
+  int ieform;
+  int lun;
+
+  int n;
+  int npar;
+  int ng;
+
+  Array<int> info;
+  Array<int> iwork;
+  Array<int> ipar;
+  Array<int> jroot;
+
+  Array<double> rwork;
+  Array<double> rpar;
+
+  Matrix y;
+  Matrix ydot;
+
+  double abs_tol;
+  double rel_tol;
+
+  double *py;
+  double *pydot;
+  int *pinfo;
+  int *piwork;
+  double *prwork;
+  double *prpar;
+  int *pipar;
+  int *pjroot;
+
+  void init_work_size (int);
+
+  void integrate (double t);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

--- a/liboctave/DASSL.cc
+++ b/liboctave/DASSL.cc
@@ -181,7 +181,7 @@
 
   // Fix up the matrix of partial derivatives for dassl.
 
-  tmp_dfdx = tmp_dfdx +  cj * tmp_dfdxdot;
+  tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot;
 
   for (int j = 0; j < nn; j++)
     for (int i = 0; i < nn; i++)

--- a/liboctave/Makefile.in
+++ b/liboctave/Makefile.in
@@ -44,18 +44,19 @@
 	vx-rv-cs.h vx-s-ccv.h vx-s-crv.h \
 	vx-rv-crv.h vx-cv-ccv.h vx-crv-rv.h vx-ccv-cv.h
 
-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 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 \
-	lo-utils.h mach-info.h oct-alloc.h oct-cmplx.h oct-env.h \
-	oct-fftw.h oct-getopt.h oct-group.h oct-kpse.h oct-passwd.h \
-	oct-rl-edit.h oct-rl-hist.h oct-shlib.h oct-syscalls.h oct-time.h \
-	pathlen.h pathsearch.h prog-args.h statdefs.h str-vec.h sun-utils.h \
-	sysdir.h systime.h syswait.h \
+INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.h DAERT.h DAERTFunc.h \
+	DASPK.h DASRT.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 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 lo-utils.h mach-info.h \
+	oct-alloc.h oct-cmplx.h oct-env.h oct-fftw.h oct-getopt.h \
+	oct-group.h oct-kpse.h oct-passwd.h oct-rl-edit.h \
+	oct-rl-hist.h oct-shlib.h oct-syscalls.h oct-time.h \
+	pathlen.h pathsearch.h prog-args.h statdefs.h str-vec.h\
+	sun-utils.h sysdir.h systime.h syswait.h \
 	$(MATRIX_INC) \
 	$(MX_OP_INC) \
 	$(VX_OP_INC)
@@ -85,10 +86,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 ODES.cc ODESSA.cc Quad.cc \
-	Range.cc data-conv.cc dir-ops.cc \
+LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc \
+	DASPK.cc DASRT.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 \

new file mode 100644
--- /dev/null
+++ b/liboctave/base-dae.h
@@ -0,0 +1,86 @@
+/*
+
+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_base_dae_h)
+#define octave_base_dae_h 1
+
+#include "base-de.h"
+
+class
+base_diff_alg_eqn : public base_diff_eqn
+{
+public:
+
+  base_diff_alg_eqn (void)
+    : base_diff_eqn (), xdot () { }
+
+  base_diff_alg_eqn (const ColumnVector& xx, double tt)
+    : base_diff_eqn (xx, tt), xdot (xx.length (), 0.0) { }
+
+  base_diff_alg_eqn (const ColumnVector& xx, const ColumnVector& xxdot,
+		     double tt)
+    : base_diff_eqn (xx, tt), xdot (xxdot) { }
+
+  base_diff_alg_eqn (const base_diff_alg_eqn& a)
+    : base_diff_eqn (a), xdot (a.xdot) { }
+
+  virtual ~base_diff_alg_eqn (void) { }
+
+  base_diff_alg_eqn& operator = (const base_diff_alg_eqn& a)
+    {
+      if (this != &a)
+	{
+	  base_diff_eqn::operator = (a);
+	  xdot = a.xdot;
+	}
+      return *this;
+    }
+
+  void initialize (const ColumnVector& x0, double t0)
+    {
+      base_diff_eqn::initialize (x0, t0);
+      xdot.resize (x0.length (), 0.0);
+      force_restart ();
+    }
+
+  void initialize (const ColumnVector& x0, const ColumnVector& xdot0,
+		   double t0)
+    {
+      base_diff_eqn::initialize (x0, t0);
+      xdot = xdot0;
+      force_restart ();
+    }
+
+  ColumnVector state_derivative (void) { return xdot; }
+
+protected:
+
+  ColumnVector xdot;
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@
+2002-07-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* DLD-FUNCTIONS/dasrt.cc: New file.
+	* Makefile.in (DLD_XSRC): Add it to the list.
+
 2002-07-12  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* lex.l (@): Handle new token.

new file mode 100644
--- /dev/null
+++ b/src/DLD-FUNCTIONS/dasrt.cc
@@ -0,0 +1,788 @@
+/*
+
+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 <iostream.h>
+
+#include "DASRT.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 "parse.h"
+#include "unwind-prot.h"
+#include "utils.h"
+#include "variables.h"
+
+// Global pointers for user defined function required by dassl.
+static octave_function *dasrt_f;
+static octave_function *dasrt_j;
+static octave_function *dasrt_cf;
+
+static DASRT_options dasrt_opts;
+
+// Is this a recursive call?
+static int call_depth = 0;
+
+static ColumnVector
+dasrt_user_f (const ColumnVector& x, const ColumnVector& xprime,
+	       double t, int& ires)
+{
+  ColumnVector retval;
+
+  octave_value_list args;
+
+  int n = x.length ();
+
+  if (n > 1)
+    {
+      args(0) = x;
+      args(1) = xprime;
+    }
+  else if (n == 1)
+    {
+      args(0) = x(0);
+      args(1) = xprime(0);
+    }
+  else
+    {
+      args(0) = Matrix ();
+      args(1) = Matrix ();
+    }
+
+  args(2) = t;
+
+  if (dasrt_f)
+    {
+      octave_value_list tmp = dasrt_f->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("dasrt");
+	  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 ("dasrt");
+	}
+      else
+	gripe_user_supplied_eval ("dasrt");
+    }
+
+  return retval;
+}
+
+static ColumnVector
+dasrt_user_cf (const ColumnVector& x, double t)
+{
+  ColumnVector retval;
+
+  octave_value_list args;
+
+  int n = x.length ();
+
+  if (n > 1)
+    args(0) = x;
+  else if (n == 1)
+    args(0) = x(0);
+  else
+    args(0) = Matrix ();
+
+  args(1) = t;
+
+  if (dasrt_cf)
+    {
+      octave_value_list tmp = dasrt_cf->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("dasrt");
+	  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 ("dasrt");
+	}
+      else
+	gripe_user_supplied_eval ("dasrt");
+    }
+
+  return retval;
+}
+
+static ColumnVector
+dasrt_dumb_cf (const ColumnVector& x, double t)
+{
+  ColumnVector retval (1, 1.0);
+  return retval;
+}
+
+#if 0
+static Matrix
+dasrt_user_mf (double t, const ColumnVector& x, const ColumnVector& xprime,
+		const double& cj, octave_function *mf)
+{
+  Matrix retval;
+
+  if (mf)
+    {
+      octave_value_list args;
+
+      int n = x.length ();
+
+      if (n > 1)
+        {
+	  args(0) = x;
+	  args(1) = xprime;
+	  args(3) = cj;
+        }
+      else if (n == 1)
+        {
+	  args(0) = x(0);
+	  args(1) = xprime(0);
+	  args(3) = cj;
+        }
+      else
+        {
+	  args(0) = Matrix ();
+	  args(1) = Matrix ();
+	  args(3) = Matrix ();
+        }
+
+      args(2) = t;
+
+      octave_value_list tmp = mf->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("dasrt");
+	  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 ("dasrt");
+	}
+      else
+	gripe_user_supplied_eval ("dasrt");
+    }
+
+  return retval;
+}
+
+static Matrix
+dasrt_user_j (const ColumnVector& x, const ColumnVector& xprime, double t)
+{
+  return dasrt_user_mf (t, x, xprime, dasrt_j);
+}
+#endif
+
+#define DASRT_ABORT \
+  do \
+    { \
+      unwind_protect::run_frame ("Fdasrt"); \
+      return retval; \
+    } \
+  while (0)
+
+#define DASRT_ABORT1(msg) \
+  do \
+    { \
+      ::error ("dasrt: " ## msg); \
+      DASRT_ABORT; \
+    } \
+  while (0)
+
+#define DASRT_ABORT2(fmt, arg) \
+  do \
+    { \
+      ::error ("dasrt: " ## fmt, arg); \
+      DASRT_ABORT; \
+    } \
+  while (0)
+
+DEFUN_DLD (dasrt, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{x}, @var{xdot}, @var{t}] =} dasrt (@var{fj} [, @var{g}], @var{x_0}, @var{xdot_0}, @var{t_out} [, @var{t_crit}])\n\
+Solve a system of differential/algebraic equations with functional\n\
+stopping criteria.\n\
+\n\
+The function to be integrated must be of the form:\n\
+@example\n\
+@var{res} = f (@var{x}, @var{xdot}, @var{t}) = 0\n\
+@end example\n\
+\n\
+The stopping condition must be of the form:\n\
+\n\
+@example\n\
+@var{res} = g (@var{x}, @var{t}) = 0\n\
+@end example\n\
+\n\
+The Jacobian (if included) must be of the form:\n\
+\n\
+@example\n\
+@var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{cj})\n\
+   =  df/dx + cj*df/dxdot\n\
+@end example\n\
+\n\
+@noindent\n\
+The following inputs are entered:\n\
+\n\
+@table @var\n\
+@item fj\n\
+The string vector containing either @var{f} or both @var{f} and @var{j}.\n\
+\n\
+@item f\n\
+The function to be integrated.\n\
+\n\
+@item g\n\
+The function with the stopping conditions.\n\
+\n\
+@item j\n\
+The optional Jacobian function.  If not included, it will be approximated\n\
+by finite differences.\n\
+\n\
+@item x_0\n\
+The initial state.\n\
+\n\
+@item xdot_0\n\
+The time derivative of the initial state.\n\
+\n\
+@item t_out\n\
+The times at which the solution should be returned.  This vector should\n\
+include the initial time a the first value.\n\
+\n\
+@item t_crit\n\
+The times at which the function is non-smooth or poorly behaved.\n\
+\n\
+@end table\n\
+\n\
+@noindent\n\
+The following outputs are returned:\n\
+\n\
+@table @var\n\
+@item x\n\
+The states at each time in @var{t}.\n\
+\n\
+@item xdot\n\
+The time derivative of the states at each time in @var{t}.\n\
+\n\
+@item t\n\
+All the times in the requested output time vector up to the stopping\n\
+criteria.  The time at which the stopping criteria is achieved is returned\n\
+as the last time in the vector.\n\
+@end table\n\
+\n\
+You can use the function @code{dasrt_options} to set optional\n\
+parameters for @code{dasrt}.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  unwind_protect::begin_frame ("Fdasrt");
+
+  unwind_protect_int (call_depth);
+  call_depth++;
+
+  if (call_depth > 1)
+    DASRT_ABORT1 ("invalid recursive call");
+
+  int argp = 0;
+  int ng = 0;
+
+  int nargin = args.length ();
+
+  if (nargin < 4 || nargin > 6)
+    {
+      print_usage ("dasrt");
+      unwind_protect::run_frame ("Fdasrt");
+      return retval;
+    }
+
+  dasrt_f = 0;
+  dasrt_j = 0;
+  dasrt_cf = 0;
+
+  // Check all the arguments.  Are they the right animals?
+
+  // Here's where I take care of f and j in one shot:
+
+  octave_value f_arg = args(0);
+
+  switch (f_arg.rows ())
+    {
+    case 1:
+      dasrt_f = extract_function
+	(args(0), "dasrt", "__dasrt_fcn__",
+	 "function res = __dasrt_fcn__ (x, xdot, t) res = ",
+	 "; endfunction");
+      break;
+      
+    case 2:
+      {
+	string_vector tmp = args(0).all_strings ();
+	
+	if (! error_state)
+	  {
+	    dasrt_f = extract_function
+	      (tmp(0), "dasrt", "__dasrt_fcn__",
+	       "function res = __dasrt_fcn__ (x, xdot, t) res = ",
+	       "; endfunction");
+	    
+	    if (dasrt_f)
+	      {
+		dasrt_j = extract_function
+		  (tmp(1), "dasrt", "__dasrt_jac__",
+		   "function jac = __dasrt_jac__ (x, xdot, t, cj) jac = ",
+		   "; endfunction");
+		
+		if (! dasrt_j)
+		  dasrt_f = 0;
+	      }
+	  }
+      }
+      break;
+      
+    default:
+      DASRT_ABORT1
+	("first arg should be a string or 2-element string array");
+    }
+  
+  if (error_state || (! dasrt_f))
+    DASRT_ABORT;
+  
+  DAERTFunc func (dasrt_user_f);
+  
+  argp++;
+  double t0;
+  ColumnVector x0;
+  
+  if (args(1).is_string ())
+    {
+      dasrt_cf = is_valid_function (args(1), "dasrt", true);
+
+      if (! dasrt_cf)
+	  DASRT_ABORT1 ("expecting function name as argument 2");
+      else
+        {
+	  octave_value_list blah;
+	  octave_value_list inputs;
+	  ColumnVector the_ts;
+	  x0 = ColumnVector (args(2).vector_value ());
+
+	  if (error_state)
+	    DASRT_ABORT1 ("bad x0");
+
+	  the_ts = ColumnVector (args(4).vector_value ());
+
+	  if (error_state)
+	    DASRT_ABORT1 ("bad tout");
+
+	  t0 = the_ts(0);
+
+	  inputs(0) = x0;
+	  inputs(1) = t0;
+
+	  string g_name = args(1).string_value();
+
+	  if (error_state)
+	    DASRT_ABORT1 ("bad function name of stopping condition");
+
+	  blah = feval (g_name, inputs, 2);
+	  ColumnVector gout = ColumnVector (blah(0).vector_value ());
+
+	  ng = gout.length();
+	  int testg;
+	  for (testg = 0; testg < ng; testg++)
+	    {
+	      if (gout(testg) == 0)
+		DASRT_ABORT1 ("stopping condition satisfied at initial point");
+	    }
+        }
+      argp++;
+
+      func.set_constraint_function (dasrt_user_cf);
+    }
+  else
+    {
+      // Now this second argument is not a string.  It has to be x0.
+      // we call the dummy g function and set ng = 1;
+
+      x0 = ColumnVector (args(1).vector_value ());
+
+      if (error_state)
+	DASRT_ABORT1 ("bad x0");
+
+      func.set_constraint_function (dasrt_dumb_cf);
+
+      ng = 1;
+    }
+
+  ColumnVector state (args(argp++).vector_value ());
+
+  if (error_state)
+    DASRT_ABORT2 ("expecting state vector as argument %d", argp);
+
+  ColumnVector stateprime (args(argp++).vector_value());
+
+  if (error_state)
+    DASRT_ABORT2 
+       ("expecting time derivative of state vector as argument %d", argp);
+
+  ColumnVector old_out_times (args(argp++).vector_value ());
+
+  if (error_state)
+    DASRT_ABORT2
+	("expecting output time vector as %s argument %d", argp);
+
+  double tzero = old_out_times (0);
+
+  int ol = old_out_times.length ();
+  int ijk;
+
+  ColumnVector out_times (ol-1, 0.0);
+
+  for (ijk = 1; ijk < ol; ijk++)
+    out_times(ijk-1) = old_out_times(ijk);
+
+  ColumnVector crit_times;
+
+  bool crit_times_set = false;
+
+  if (argp < nargin)
+    {
+      crit_times = ColumnVector (args(argp++).vector_value ());
+
+      if (error_state)
+	DASRT_ABORT2
+	  ("expecting critical time vector as argument %d", argp);
+
+      crit_times_set = true;
+    }
+
+#if 0
+  if (dasrt_j)
+    func.set_jacobian_function (dasrt_user_j);
+#endif
+
+  DASRT_result output;
+
+  DASRT dae = DASRT (ng, state, stateprime, tzero, func);
+
+  dae.copy (dasrt_opts);
+  dae.set_ng (ng);
+
+  if (error_state)
+    DASRT_ABORT1 ("something is wrong");
+
+  if (crit_times_set)
+    output = dae.integrate (out_times, crit_times);
+  else
+    output = dae.integrate (out_times);
+
+  if (! error_state)
+    {
+      ColumnVector old_output_times = output.times ();
+
+      Matrix old_output_deriv = output.deriv ();
+      Matrix old_output_state = output.state ();
+
+      int lstuff = old_output_times.length ();
+      int lstate = x0.length ();
+
+      ColumnVector output_times (lstuff+1, 0.0);
+
+      Matrix output_deriv (lstuff+1, lstate, 0.0);
+      Matrix output_state (lstuff+1, lstate, 0.0);
+
+      output_times(0) = tzero;
+
+      for (ijk = 0; ijk < lstate; ijk++)
+	{
+	  output_deriv(0,ijk) = stateprime(ijk);
+	  output_state(0,ijk) = state(ijk);
+	}
+
+      for (ijk = 0; ijk < lstuff; ijk++)
+	{
+	  output_times(ijk+1) = old_output_times(ijk);
+
+	  for (int lmnop=0; lmnop < lstate; lmnop++)
+	    {
+	      output_deriv(ijk+1,lmnop) = old_output_deriv(ijk,lmnop);
+	      output_state(ijk+1,lmnop) = old_output_state(ijk,lmnop);
+	    }
+	}
+
+      retval(2) = output_times;
+      retval(1) = output_deriv;
+      retval(0) = output_state;
+    }
+  else
+    {
+      DASRT_ABORT1("something wicked has occurred!");
+      // print_usage ("dasrt");
+    }
+
+  unwind_protect::run_frame ("Fdasrt");
+
+  return retval;
+}
+
+typedef void (DASRT_options::*d_set_opt_mf) (double);
+typedef void (DASRT_options::*i_set_opt_mf) (int);
+typedef double (DASRT_options::*d_get_opt_mf) (void);
+typedef int (DASRT_options::*i_get_opt_mf) (void);
+
+#define MAX_TOKENS 3
+
+struct DASRT_OPTIONS
+{
+  const char *keyword;
+  const char *kw_tok[MAX_TOKENS + 1];
+  int min_len[MAX_TOKENS + 1];
+  int min_toks_to_match;
+  d_set_opt_mf d_set_fcn;
+  i_set_opt_mf i_set_fcn;
+  d_get_opt_mf d_get_fcn;
+  i_get_opt_mf i_get_fcn;
+};
+
+static DASRT_OPTIONS dasrt_option_table [] =
+{
+  { "absolute tolerance",
+    { "absolute", "tolerance", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    &DASRT_options::set_absolute_tolerance, 0,
+    &DASRT_options::absolute_tolerance, 0, },
+
+  { "initial step size",
+    { "initial", "step", "size", 0, },
+    { 1, 0, 0, 0, }, 1,
+    &DASRT_options::set_initial_step_size, 0,
+    &DASRT_options::initial_step_size, 0, },
+
+  { "maximum step size",
+    { "maximum", "step", "size", 0, },
+    { 2, 0, 0, 0, }, 1,
+    &DASRT_options::set_maximum_step_size, 0,
+    &DASRT_options::maximum_step_size, 0, },
+
+  { "minimum step size",
+    { "minimum", "step", "size", 0, },
+    { 2, 0, 0, 0, }, 1,
+    &DASRT_options::set_minimum_step_size, 0,
+    &DASRT_options::minimum_step_size, 0, },
+
+  { "relative tolerance",
+    { "relative", "tolerance", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    &DASRT_options::set_relative_tolerance, 0,
+    &DASRT_options::relative_tolerance, 0, },
+
+  { "step limit",
+    { "step", "limit", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    0, &DASRT_options::set_step_limit,
+    0, &DASRT_options::step_limit, },
+
+  { 0,
+    { 0, 0, 0, 0, },
+    { 0, 0, 0, 0, }, 0,
+    0, 0, 0, 0, },
+};
+
+static void
+print_dasrt_option_list (ostream& os)
+{
+  print_usage ("dasrt_options", 1);
+
+  os << "\n"
+     << "Options for dasrt include:\n\n"
+     << "  keyword                                  value\n"
+     << "  -------                                  -----\n\n";
+
+  DASRT_OPTIONS *list = dasrt_option_table;
+
+  const char *keyword;
+  while ((keyword = list->keyword) != 0)
+    {
+      os.form ("  %-40s ", keyword);
+
+      if (list->d_get_fcn)
+	{
+	  double val = (dasrt_opts.*list->d_get_fcn) ();
+	  if (val < 0.0)
+	    os << "computed automatically";
+	  else
+	    os << val;
+	}
+      else
+	{
+	  int val = (dasrt_opts.*list->i_get_fcn) ();
+	  if (val < 0)
+	    os << "computed automatically";
+	  else
+	    os << val;
+	}
+      os << "\n";
+      list++;
+    }
+
+  os << "\n";
+}
+
+static void
+set_dasrt_option (const string& keyword, double val)
+{
+  DASRT_OPTIONS *list = dasrt_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->d_set_fcn)
+	    (dasrt_opts.*list->d_set_fcn) (val);
+	  else
+	    {
+	      if (xisnan (val))
+		{
+		  error ("dasrt_options: %s: expecting integer, found NaN",
+			 keyword.c_str ());
+		}
+	      else
+		(dasrt_opts.*list->i_set_fcn) (NINT (val));
+	    }
+	  return;
+	}
+      list++;
+    }
+
+  warning ("dasrt_options: no match for `%s'", keyword.c_str ());
+}
+
+static octave_value_list
+show_dasrt_option (const string& keyword)
+{
+  octave_value retval;
+
+  DASRT_OPTIONS *list = dasrt_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->d_get_fcn)
+	    {
+	      double val = (dasrt_opts.*list->d_get_fcn) ();
+	      if (val < 0.0)
+		retval = "computed automatically";
+	      else
+		retval = val;
+	    }
+	  else
+	    {
+	      int val = (dasrt_opts.*list->i_get_fcn) ();
+	      if (val < 0)
+		retval = "computed automatically";
+	      else
+		retval = static_cast<double> (val);
+	    }
+
+	  return retval;
+	}
+      list++;
+    }
+
+  warning ("dasrt_options: no match for `%s'", keyword.c_str ());
+
+  return retval;
+}
+
+DEFUN_DLD (dasrt_options, args, ,
+  "dasrt_options (KEYWORD, VALUE)\n\
+\n\
+Set or show options for dasrt.  Keywords may be abbreviated\n\
+to the shortest match.")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+
+  if (nargin == 0)
+    {
+      print_dasrt_option_list (octave_stdout);
+      return retval;
+    }
+  else if (nargin == 1 || nargin == 2)
+    {
+      string keyword = args(0).string_value ();
+
+      if (! error_state)
+	{
+	  if (nargin == 1)
+	    return show_dasrt_option (keyword);
+	  else
+	    {
+	      double val = args(1).double_value ();
+
+	      if (! error_state)
+		{
+		  set_dasrt_option (keyword, val);
+		  return retval;
+		}
+	    }
+	}
+    }
+
+  print_usage ("dasrt_options");
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/

--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -40,7 +40,7 @@
 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 \
+	daspk.cc dasrt.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 \