Mercurial > hg > octave-nkf

changeset 4049:a35a3c5d4740

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-08-16 08:54:31 by jwe]
author jwe
date Fri, 16 Aug 2002 08:54:31 +0000
parents c9991c59cf5c
children 6481f41a79f3
files liboctave/ChangeLog liboctave/DAEFunc.h liboctave/DAERTFunc.h liboctave/DASPK-opts.in liboctave/DASPK.cc liboctave/DASPK.h liboctave/DASRT-opts.in liboctave/DASRT.cc liboctave/DASRT.h liboctave/DASSL-opts.in liboctave/DASSL.cc liboctave/DASSL.h liboctave/LSODE-opts.in liboctave/LSODE.cc liboctave/LSODE.h liboctave/ODEFunc.h liboctave/ODESSA-opts.in liboctave/ODESSA.h liboctave/base-de.h mk-opts.pl
diffstat 20 files changed, 766 insertions(+), 747 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,5 +1,29 @@
+2002-08-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* LSODE.h (rel_tol, abs_tol, px, pabs_tol, piwork, prwork, itol):
+	New data members.
+	(LSODE::sanity_checked): Delete unused data member.
+
+	* DASPKL.h (initialized, abs_tol, rel_tol, px, pxdot, pabs_tol,
+	prel_tol, pinfo, piwork, prwork): New data members.
+	* DASSL.h (DASSL): Likewise.
+
+	* DASRT.h (DASRT::sanity_checked): Delete unused data member.
+
+	* DASRT.cc (DASRT::integrate (double)): Better handling of
+	initialization, changes in options, etc.
+	* DASPK.cc (DASPK::do_integrate): Likewise.
+	* DASSL.cc (DASSL::do_integrate): Likewise.
+	* LSODE.cc (LSODE::do_integrate): Likewise.
+
 2002-08-15  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* DAEFunc.h (DAEFunc::reset): New data member.
+	* DAERTFunc.h (DAERTFunc::reset): Likewise.
+
+	* base-de.h (base_diff_eqn::set_stop_time): Force restart here.
+	(base_diff_eqn::clear_stop_time): Likewise.
+
 	* DASSL.cc (DASSL::do_integrate (double)): Handle more optoins.
 	* DASPK.cc (DASPK::do_integrate (double)): Likewise.
 
--- a/liboctave/DAEFunc.h
+++ b/liboctave/DAEFunc.h
@@ -44,16 +44,16 @@
 				double t, double cj);
 
   DAEFunc (void)
-    : fun (0), jac (0) { }
+    : fun (0), jac (0), reset (true) { }
 
   DAEFunc (DAERHSFunc f)
-    : fun (f), jac (0) { }
+    : fun (f), jac (0), reset (true) { }
 
   DAEFunc (DAERHSFunc f, DAEJacFunc j)
-    : fun (f), jac (j) { }
+    : fun (f), jac (j), reset (true) { }
 
   DAEFunc (const DAEFunc& a)
-    : fun (a.fun), jac (a.jac) { }
+    : fun (a.fun), jac (a.jac), reset (a.reset) { }
 
   DAEFunc& operator = (const DAEFunc& a)
     {
@@ -61,6 +61,7 @@
 	{
 	  fun = a.fun;
 	  jac = a.jac;
+	  reset = a.reset;
 	}
       return *this;
     }
@@ -72,6 +73,7 @@
   DAEFunc& set_function (DAERHSFunc f)
     {
       fun = f;
+      reset = true;
       return *this;
     }
 
@@ -80,6 +82,7 @@
   DAEFunc& set_jacobian_function (DAEJacFunc j)
     {
       jac = j;
+      reset = true;
       return *this;
     }
 
@@ -87,6 +90,13 @@
 
   DAERHSFunc fun;
   DAEJacFunc jac;
+
+  // This variable is TRUE when this object is constructed, and also
+  // after any internal data has changed.  Derived classes may use
+  // this information (and change it) to know when to (re)initialize
+  // their own internal data related to this object.
+
+  bool reset;
 };
 
 #endif
--- a/liboctave/DAERTFunc.h
+++ b/liboctave/DAERTFunc.h
@@ -33,22 +33,22 @@
   typedef ColumnVector (*DAERTConstrFunc) (const ColumnVector& x, double t);
 
   DAERTFunc (void)
-    : DAEFunc (), constr (0) { }
+    : DAEFunc (), constr (0), reset (true) { }
 
   DAERTFunc (DAERHSFunc f)
-    : DAEFunc (f), constr (0) { }
+    : DAEFunc (f), constr (0), reset (true) { }
 
   DAERTFunc (DAERHSFunc f, DAEJacFunc j)
-    : DAEFunc (f, j), constr (0) { }
+    : DAEFunc (f, j), constr (0), reset (true) { }
 
   DAERTFunc (DAERHSFunc f, DAERTConstrFunc cf)
-    : DAEFunc (f), constr (cf) { }
+    : DAEFunc (f), constr (cf), reset (true) { }
 
   DAERTFunc (DAERHSFunc f, DAERTConstrFunc cf, DAEJacFunc j)
-    : DAEFunc (f, j), constr (cf) { }
+    : DAEFunc (f, j), constr (cf), reset (true) { }
 
   DAERTFunc (const DAERTFunc& a)
-    : DAEFunc (a), constr (a.constr) { }
+    : DAEFunc (a), constr (a.constr), reset (a.reset) { }
 
   DAERTFunc& operator = (const DAERTFunc& a)
     {
@@ -56,6 +56,7 @@
 	{
 	  DAEFunc::operator = (a);
 	  constr = a.constr;
+	  reset = a.reset;
 	}
       return *this;
     }
@@ -67,12 +68,20 @@
   DAERTFunc& set_constraint_function (DAERTConstrFunc cf)
     {
       constr = cf;
+      reset = true;
       return *this;
     }
 
 protected:
 
   DAERTConstrFunc constr;
+
+  // This variable is TRUE when this object is constructed, and also
+  // after any internal data has changed.  Derived classes may use
+  // this information (and change it) to know when to (re)initialize
+  // their own internal data related to this object.
+
+  bool reset;
 };
 
 #endif
--- a/liboctave/DASPK-opts.in
+++ b/liboctave/DASPK-opts.in
@@ -15,10 +15,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
@@ -35,10 +36,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
@@ -62,10 +64,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = val;
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
@@ -89,10 +92,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = val;
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
@@ -153,4 +157,3 @@
   INIT_VALUE = "0"
   SET_EXPR = "val"
 END_OPTION
-
--- a/liboctave/DASPK.cc
+++ b/liboctave/DASPK.cc
@@ -64,47 +64,7 @@
 static DAEFunc::DAEJacFunc user_jac;
 static int nn;
 
-DASPK::DASPK (void) : DAE ()
-{
-  sanity_checked = 0;
-
-  info.resize (20);
-
-  for (int i = 0; i < 20; i++)
-    info.elem (i) = 0;
-}
-
-DASPK::DASPK (const ColumnVector& state, double time, DAEFunc& f)
-  : DAE (state, time, f)
-{
-  n = size ();
-
-  sanity_checked = 0;
-
-  info.resize (20);
-
-  for (int i = 0; i < 20; i++)
-    info.elem (i) = 0;
-}
-
-DASPK::DASPK (const ColumnVector& state, const ColumnVector& deriv,
-	  double time, DAEFunc& f)
-  : DAE (state, deriv, time, f)
-{
-  n = size ();
-
-  DAEFunc::set_function (f.function ());
-  DAEFunc::set_jacobian_function (f.jacobian_function ());
-
-  sanity_checked = 0;
-
-  info.resize (20);
-
-  for (int i = 0; i < 20; i++)
-    info.elem (i) = 0;
-}
-
-int
+static int
 ddaspk_f (const double& time, const double *state, const double *deriv,
 	  const double&, double *delta, int& ires, double *, int *)
 {
@@ -137,7 +97,7 @@
 //NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT,
 //C                          WP, IWP, B, EPLIN, IER, RPAR, IPAR)
 
-int
+static int
 ddaspk_psol (const int& neq, const double& time, const double *state,
 	     const double *deriv, const double *savr,
 	     const double& cj, const double *wght, double *wp,
@@ -149,7 +109,7 @@
 }
 
 
-int
+static int
 ddaspk_j (const double& time, const double *state, const double *deriv,
 	  double *pd, const double& cj, double *, int *)
 {
@@ -181,178 +141,220 @@
 
   ColumnVector retval;
 
-  if (restart)
+  if (! initialized || restart || DAEFunc::reset|| DASPK_options::reset)
     {
-      restart = false;
-      info.elem (0) = 0;
-    }
+      integration_error = false;
+
+      initialized = true;
+
+      info.resize (20);
 
-  liw = 40 + n;
-  if (info(9) == 1 || info(9) == 3)
-    liw += n;
-  if (info (10) == 1 || info(15) == 1)
-    liw += n;
+      for (int i = 0; i < 20; i++)
+	info(i) = 0;
+
+      pinfo = info.fortran_vec ();
 
-  lrw = 50 + 9*n;
-  if (info(5) == 0)
-    lrw += n*n;
-  if (info(15) == 1)
-    lrw += n;
+      int n = size ();
 
-  if (iwork.length () != liw)
-    iwork.resize (liw);
+      nn = n;
+
+      info(0) = 0;
 
-  if (rwork.length () != lrw)
-    rwork.resize (lrw);
-
-  integration_error = false;
+      if (stop_time_set)
+	{
+	  rwork(0) = stop_time;
+	  info(3) = 1;
+	}
+      else
+	info(3) = 0;
 
-  if (DAEFunc::jacobian_function ())
-    info.elem (4) = 1;
-  else
-    info.elem (4) = 0;
+      px = x.fortran_vec ();
+      pxdot = xdot.fortran_vec ();
+
+      // DAEFunc
+
+      user_fun = DAEFunc::function ();
+      user_jac = DAEFunc::jacobian_function ();
 
-  double *px    = x.fortran_vec ();
-  double *pxdot = xdot.fortran_vec ();
+      if (user_fun)
+	{
+	  int ires = 0;
 
-  nn = n;
-  user_fun = DAEFunc::fun;
-  user_jac = DAEFunc::jac;
+	  ColumnVector res = (*user_fun) (x, xdot, t, ires);
 
-  if (! sanity_checked)
-    {
-      int ires = 0;
+	  if (res.length () != x.length ())
+	    {
+	      (*current_liboctave_error_handler)
+		("daspk: inconsistent sizes for state and residual vectors");
 
-      ColumnVector res = (*user_fun) (x, xdot, t, ires);
-
-      if (res.length () != x.length ())
+	      integration_error = true;
+	      return retval;
+	    }
+	}
+      else
 	{
 	  (*current_liboctave_error_handler)
-	    ("daspk: inconsistent sizes for state and residual vectors");
+	    ("daspk: no user supplied RHS subroutine!");
+
+	  integration_error = true;
+	  return retval;
+	}
+  
+      info(4) = user_jac ? 1 : 0;
+
+      DAEFunc::reset = false;
+
+      // DASPK_options
+
+      Array<double> abs_tol = absolute_tolerance ();
+      Array<double> rel_tol = relative_tolerance ();
+
+      int abs_tol_len = abs_tol.length ();
+      int rel_tol_len = rel_tol.length ();
+
+      if (abs_tol_len == 1 && rel_tol_len == 1)
+	{
+	  info(1) = 0;
+	}
+      else if (abs_tol_len == n && rel_tol_len == n)
+	{
+	  info(1) = 1;
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("daspk: inconsistent sizes for tolerance arrays");
 
 	  integration_error = true;
 	  return retval;
 	}
 
-      sanity_checked = 1;
-    }
-  
-  if (stop_time_set)
-    {
-      rwork.elem (0) = stop_time;
-      info.elem (3) = 1;
-    }
-  else
-    info.elem (3) = 0;
-
-  Array<double> abs_tol = absolute_tolerance ();
-  Array<double> rel_tol = relative_tolerance ();
-
-  int abs_tol_len = abs_tol.length ();
-  int rel_tol_len = rel_tol.length ();
-
-  if (abs_tol_len == 1 && rel_tol_len == 1)
-    {
-      info.elem (1) = 0;
-    }
-  else if (abs_tol_len == n && rel_tol_len == n)
-    {
-      info.elem (1) = 1;
-    }
-  else
-    {
-      (*current_liboctave_error_handler)
-	("daspk: inconsistent sizes for tolerance arrays");
+      double hmax = maximum_step_size ();
+      if (hmax >= 0.0)
+	{
+	  rwork(1) = hmax;
+	  info(6) = 1;
+	}
+      else
+	info(6) = 0;
 
-      integration_error = true;
-      return retval;
-    }
-
-  double hmax = maximum_step_size ();
-  if (hmax >= 0.0)
-    {
-      rwork.elem (1) = hmax;
-      info.elem (6) = 1;
-    }
-  else
-    info.elem (6) = 0;
-
-  double h0 = initial_step_size ();
-  if (h0 >= 0.0)
-    {
-      rwork.elem (2) = h0;
-      info.elem (7) = 1;
-    }
-  else
-    info.elem (7) = 0;
-
-  int maxord = maximum_order ();
-  if (maxord >= 0)
-    {
-      if (maxord > 0 && maxord < 6)
+      double h0 = initial_step_size ();
+      if (h0 >= 0.0)
 	{
-	  info(8) = 1;
-	  iwork(2) = maxord;
+	  rwork(2) = h0;
+	  info(7) = 1;
 	}
       else
+	info(7) = 0;
+
+      int maxord = maximum_order ();
+      if (maxord >= 0)
 	{
+	  if (maxord > 0 && maxord < 6)
+	    {
+	      info(8) = 1;
+	      iwork(2) = maxord;
+	    }
+	  else
+	    {
+	      (*current_liboctave_error_handler)
+		("daspk: invalid value for maximum order");
+	      integration_error = true;
+	      return retval;
+	    }
+	}
+
+      int eiq = enforce_inequality_constraints ();
+      switch (eiq)
+	{
+	case 1:
+	case 3:
+	  {
+	    Array<int> ict = inequality_constraint_types ();
+
+	    if (ict.length () == n)
+	      {
+		for (int i = 0; i < n; i++)
+		  {
+		    int val = ict(i);
+		    if (val < -2 || val > 2)
+		      {
+			(*current_liboctave_error_handler)
+			  ("daspk: invalid value for inequality constraint type");
+			integration_error = true;
+			return retval;
+		      }
+		    iwork(40+i) = val;
+		  }
+	      }
+	    else
+	      {
+		(*current_liboctave_error_handler)
+		  ("daspk: inequality constraint types size mismatch");
+		integration_error = true;
+		return retval;
+	      }
+	  }
+	  // Fall through...
+
+	case 2:
+	  info(9) = eiq;
+	  break;
+
+	default:
 	  (*current_liboctave_error_handler)
-	    ("daspk: invalid value for maximum order");
+	    ("daspk: invalid value for enforce inequality constraints option");
 	  integration_error = true;
 	  return retval;
 	}
-    }
+
+      int ccic = compute_consistent_initial_condition ();
+      if (ccic)
+	{
+	  if (ccic == 1)
+	    {
+	      // XXX FIXME XXX -- this code is duplicated below.
+
+	      Array<int> av = algebraic_variables ();
 
-  int eiq = enforce_inequality_constraints ();
-  switch (eiq)
-    {
-    case 1:
-    case 3:
-      {
-	Array<int> ict = inequality_constraint_types ();
+	      if (av.length () == n)
+		{
+		  int lid;
+		  if (eiq == 0 || eiq == 2)
+		    lid = 40;
+		  else if (eiq == 1 || eiq == 3)
+		    lid = 40 + n;
+		  else
+		    abort ();
 
-	if (ict.length () == n)
-	  {
-	    for (int i = 0; i < n; i++)
-	      {
-		int val = ict(i);
-		if (val < -2 || val > 2)
-		  {
-		    (*current_liboctave_error_handler)
-		      ("daspk: invalid value for inequality constraint type");
-		    integration_error = true;
-		    return retval;
-		  }
-		iwork(40+i) = val;
-	      }
-	  }
-	else
-	  {
-	    (*current_liboctave_error_handler)
-	      ("daspk: inequality constraint types size mismatch");
-	    integration_error = true;
-	    return retval;
-	  }
-      }
-      // Fall through...
+		  for (int i = 0; i < n; i++)
+		    iwork(lid+i) = av(i) ? -1 : 1;
+		}
+	      else
+		{
+		  (*current_liboctave_error_handler)
+		    ("daspk: algebraic variables size mismatch");
+		  integration_error = true;
+		  return retval;
+		}
+	    }
+	  else if (ccic != 2)
+	    {
+	      (*current_liboctave_error_handler)
+		("daspk: invalid value for compute consistent initial condition option");
+	      integration_error = true;
+	      return retval;
+	    }
 
-    case 2:
-      info(9) = eiq;
-      break;
+	  info(10) = ccic;
+	}
 
-    default:
-      (*current_liboctave_error_handler)
-	("daspk: invalid value for enforce inequality constraints option");
-      integration_error = true;
-      return retval;
-    }
+      int eavfet = exclude_algebraic_variables_from_error_test ();
+      if (eavfet)
+	{
+	  info(15) = 1;
 
-  int ccic = compute_consistent_initial_condition ();
-  if (ccic)
-    {
-      if (ccic == 1)
-	{
-	  // XXX FIXME XXX -- this code is duplicated below.
+	  // XXX FIXME XXX -- this code is duplicated above.
 
 	  Array<int> av = algebraic_variables ();
 
@@ -369,102 +371,77 @@
 	      for (int i = 0; i < n; i++)
 		iwork(lid+i) = av(i) ? -1 : 1;
 	    }
+	}
+
+      if (use_initial_condition_heuristics ())
+	{
+	  Array<double> ich = initial_condition_heuristics ();
+
+	  if (ich.length () == 6)
+	    {
+	      iwork(31) = NINT (ich(0));
+	      iwork(32) = NINT (ich(1));
+	      iwork(33) = NINT (ich(2));
+	      iwork(34) = NINT (ich(3));
+
+	      rwork(13) = ich(4);
+	      rwork(14) = ich(5);
+	    }
 	  else
 	    {
 	      (*current_liboctave_error_handler)
-		("daspk: algebraic variables size mismatch");
+		("daspk: invalid initial condition heuristics option");
 	      integration_error = true;
 	      return retval;
 	    }
+
+	  info(16) = 1;
 	}
-      else if (ccic != 2)
+
+      int pici = print_initial_condition_info ();
+      switch (pici)
 	{
+	case 0:
+	case 1:
+	case 2:
+	  info(17) = pici;
+	  break;
+
+	default:
 	  (*current_liboctave_error_handler)
-	    ("daspk: invalid value for compute consistent initial condition option");
+	    ("daspk: invalid value for print initial condition info option");
 	  integration_error = true;
 	  return retval;
+	  break;
 	}
 
-      info(10) = ccic;
-    }
+      DASPK_options::reset = false;
 
-  if (exclude_algebraic_variables_from_error_test ())
-    {
-      info(15) = 1;
-
-      // XXX FIXME XXX -- this code is duplicated above.
-
-      Array<int> av = algebraic_variables ();
+      liw = 40 + n;
+      if (eiq == 1 || eiq == 3)
+	liw += n;
+      if (ccic == 1 || eavfet == 1)
+	liw += n;
 
-      if (av.length () == n)
-	{
-	  int lid;
-	  if (eiq == 0 || eiq == 2)
-	    lid = 40;
-	  else if (eiq == 1 || eiq == 3)
-	    lid = 40 + n;
-	  else
-	    abort ();
+      lrw = 50 + 9*n;
+      if (! user_jac)
+	lrw += n*n;
+      if (eavfet == 1)
+	lrw += n;
 
-	  for (int i = 0; i < n; i++)
-	    iwork(lid+i) = av(i) ? -1 : 1;
-	}
+      iwork.resize (liw);
+      rwork.resize (lrw);
+
+      piwork = iwork.fortran_vec ();
+      prwork = rwork.fortran_vec ();
+
+      restart = false;
     }
 
-  if (use_initial_condition_heuristics ())
-    {
-      Array<double> ich = initial_condition_heuristics ();
-
-      if (ich.length () == 6)
-	{
-	  iwork(31) = NINT (ich(0));
-	  iwork(32) = NINT (ich(1));
-	  iwork(33) = NINT (ich(2));
-	  iwork(34) = NINT (ich(3));
-
-	  rwork(13) = ich(4);
-	  rwork(14) = ich(5);
-	}
-      else
-	{
-	  (*current_liboctave_error_handler)
-	    ("daspk: invalid initial condition heuristics option");
-	  integration_error = true;
-	  return retval;
-	}
-
-      info(16) = 1;
-    }
+  static double *dummy = 0;
+  static int *idummy = 0;
 
-  int pici = print_initial_condition_info ();
-  switch (pici)
-    {
-    case 0:
-    case 1:
-    case 2:
-      info(17) = pici;
-      break;
-
-    default:
-      (*current_liboctave_error_handler)
-	("daspk: invalid value for print initial condition info option");
-      integration_error = true;
-      return retval;
-      break;
-    }
-
-  double *dummy = 0;
-  int *idummy = 0;
-
-  int *pinfo = info.fortran_vec ();
-  int *piwork = iwork.fortran_vec ();
-  double *prwork = rwork.fortran_vec ();
-  double *pabs_tol = abs_tol.fortran_vec ();
-  double *prel_tol = rel_tol.fortran_vec ();
-
-// again:
-
-  F77_XFCN (ddaspk, DDASPK, (ddaspk_f, n, t, px, pxdot, tout, pinfo,
+  F77_XFCN (ddaspk, DDASPK, (ddaspk_f, nn, t, px, pxdot, tout, pinfo,
 			     prel_tol, pabs_tol, istate, prwork, lrw,
 			     piwork, liw, dummy, idummy, ddaspk_j,
 			     ddaspk_psol));
@@ -546,7 +523,9 @@
 DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out)
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
@@ -589,7 +568,9 @@
 		  const ColumnVector& tcrit) 
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
--- a/liboctave/DASPK.h
+++ b/liboctave/DASPK.h
@@ -37,12 +37,14 @@
 {
 public:
 
-  DASPK (void);
+  DASPK (void) : DAE (), DASPK_options (), initialized (false) { }
 
-  DASPK (const ColumnVector& x, double time, DAEFunc& f);
+  DASPK (const ColumnVector& state, double time, DAEFunc& f)
+    : DAE (state, time, f), DASPK_options (), initialized (false) { }
 
-  DASPK (const ColumnVector& x, const ColumnVector& xdot,
-	 double time, DAEFunc& f);
+  DASPK (const ColumnVector& state, const ColumnVector& deriv,
+	 double time, DAEFunc& f)
+    : DAE (state, deriv, time, f), DASPK_options (), initialized (false) { }
 
   ~DASPK (void) { }
 
@@ -61,20 +63,26 @@
 
 private:
 
-  int n;
+  bool initialized;
+
   int liw;  
   int lrw;
-  int sanity_checked;
+
   Array<int> info;
   Array<int> iwork;
+
   Array<double> rwork;
 
-  friend int ddaspk_j (double *time, double *state, double *deriv,
-		       double *pd, double *cj, double *rpar, int *ipar);
+  Array<double> abs_tol;
+  Array<double> rel_tol;
 
-  friend int ddaspk_f (double *time, double *state, double *deriv,
-		       double *delta, int *ires, double *rpar, int *ipar);
-
+  double *px;
+  double *pxdot;
+  double *pabs_tol;
+  double *prel_tol;
+  int *pinfo;
+  int *piwork;
+  double *prwork;
 };
 
 #endif
--- a/liboctave/DASRT-opts.in
+++ b/liboctave/DASRT-opts.in
@@ -15,10 +15,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
@@ -35,10 +36,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
@@ -50,6 +52,13 @@
 END_OPTION
 
 OPTION
+  NAME = "maximum order"
+  TYPE = "int"
+  INIT_VALUE = "-1"
+  SET_EXPR = "val"
+END_OPTION
+
+OPTION
   NAME = "maximum step size"
   TYPE = "double"
   INIT_VALUE = "-1.0"
@@ -57,13 +66,6 @@
 END_OPTION
 
 OPTION
-  NAME = "minimum step size"
-  TYPE = "double"
-  INIT_VALUE = "0.0"
-  SET_EXPR = "(val >= 0.0) ? val : 0.0"
-END_OPTION
-
-OPTION
   NAME = "step limit"
   TYPE = "int"
   INIT_VALUE = "-1"
--- a/liboctave/DASRT.cc
+++ b/liboctave/DASRT.cc
@@ -48,10 +48,6 @@
 #include "f77-fcn.h"
 #include "lo-error.h"
 
-#ifndef F77_FUNC
-#define F77_FUNC(x, X) F77_FCN (x, X)
-#endif
-
 typedef int (*dasrt_fcn_ptr) (const double&, const double*, const double*,
 			      double*, int&, double*, int*);
 
@@ -142,102 +138,63 @@
   return 0;
 }
 
-DASRT::DASRT (void)
-  : DAERT ()
-{
-  initialized = false;
-
-  sanity_checked = false;
-
-  info.resize (30, 0);
-
-  ng = 0;
-
-  liw = 0;
-  lrw = 0;
-}
-
-DASRT::DASRT (const ColumnVector& state, double time, DAERTFunc& f)
-  : DAERT (state, time, f)
-{
-  n = size ();
-
-  initialized = false;
-
-  liw = 20 + n;
-  lrw = 50 + 9*n + n*n;
-
-  sanity_checked = false;
-
-  info.resize (15, 0);
-
-  DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function ();
-  
-  if (tmp_csub)
-    {
-      ColumnVector tmp = tmp_csub (state, time);
-      ng = tmp.length ();
-    }
-  else
-    ng = 0;
-
-  jroot.resize (ng, 1);
-}
-
-DASRT::DASRT (const ColumnVector& state, const ColumnVector& deriv,
-	      double time, DAERTFunc& f)
-  : DAERT (state, deriv, time, f)
-{
-  n = size ();
-
-  initialized = false;
-
-  sanity_checked = false;
-
-  info.resize (30, 0);
-
-  DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function ();
-  
-  if (tmp_csub)
-    {
-      ColumnVector tmp = tmp_csub (state, time);
-      ng = tmp.length ();
-    }
-  else
-    ng = 0;
-
-  liw = 20 + n + 1000;
-  lrw = 50 + 9*n + n*n + 1000;
-
-  jroot.resize (ng, 1);
-}
-
 void
 DASRT::integrate (double tout)
 {
   DASRT_result retval;
 
-  if (! initialized)
+  // I suppose this is the safe thing to do.  If this is the first
+  // call, or if anything about the problem has changed, we should
+  // start completely fresh.
+
+  if (! initialized || restart
+      || DAEFunc::reset || DAERTFunc::reset || DASRT_options::reset)
     {
-      info(0) = 0;
-
       integration_error = false;
 
-      user_fsub = DAEFunc::function ();
-      user_jsub = DAEFunc::jacobian_function ();
-      user_csub = DAERTFunc::constraint_function ();
+      initialized = true;
+
+      info.resize (15);
+
+      for (int i = 0; i < 15; i++)
+	info(i) = 0;
+
+      pinfo = info.fortran_vec ();
+
+      int n = size ();
+
+      nn = n;
 
-      if (user_jsub)
-	info(4) = 1;
+      liw = 20 + n;
+      lrw = 50 + 9*n + n*n;
+
+      iwork.resize (liw);
+      rwork.resize (lrw);
+
+      info(0) = 0;
+
+      if (stop_time_set)
+	{
+	  info(3) = 1;
+	  rwork(0) = stop_time;
+	}
       else
-	info(4) = 0;
+	info(3) = 0;
 
       px = x.fortran_vec ();
       pxdot = xdot.fortran_vec ();
 
-      nn = n;
+      piwork = iwork.fortran_vec ();
+      prwork = rwork.fortran_vec ();
+
+      restart = false;
 
-      if (! sanity_checked)
+      // DAEFunc
+
+      user_fsub = DAEFunc::function ();
+      user_jsub = DAEFunc::jacobian_function ();
+
+      if (user_fsub)
 	{
 	  int ires = 0;
 
@@ -246,20 +203,85 @@
 	  if (fval.length () != x.length ())
 	    {
 	      (*current_liboctave_error_handler)
-		("dassl: inconsistent sizes for state and residual vectors");
+		("dasrt: inconsistent sizes for state and residual vectors");
 
 	      integration_error = true;
 	      return;
 	    }
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("dasrt: no user supplied RHS subroutine!");
 
-	  sanity_checked = true;
+	  integration_error = true;
+	  return;
+	}
+
+      info(4) = user_jsub ? 1 : 0;
+
+      DAEFunc::reset = false;
+
+      // DAERTFunc
+
+      user_csub = DAERTFunc::constraint_function ();
+
+      if (user_csub)
+	{
+	  ColumnVector tmp = (*user_csub) (x, t);
+	  ng = tmp.length ();
 	}
-  
-      if (iwork.length () != liw)
-	iwork.resize (liw);
+      else
+	ng = 0;
+
+      jroot.resize (ng, 1);
+
+      pjroot = jroot.fortran_vec ();
+
+      DAERTFunc::reset = false;
+
+      // DASRT_options
+
+      if (maximum_step_size () >= 0.0)
+	{
+	  rwork(1) = maximum_step_size ();
+	  info(6) = 1;
+	}
+      else
+	info(6) = 0;
 
-      if (rwork.length () != lrw)
-	rwork.resize (lrw);
+      if (initial_step_size () >= 0.0)
+	{
+	  rwork(2) = initial_step_size ();
+	  info(7) = 1;
+	}
+      else
+	info(7) = 0;
+
+      int maxord = maximum_order ();
+      if (maxord >= 0)
+	{
+	  if (maxord > 0 && maxord < 6)
+	    {
+	      info(8) = 1;
+	      iwork(2) = maxord;
+	    }
+	  else
+	    {
+	      (*current_liboctave_error_handler)
+		("dassl: invalid value for maximum order");
+	      integration_error = true;
+	      return;
+	    }
+	}
+
+      if (step_limit () >= 0)
+	{
+	  info(11) = 1;
+	  iwork(18) = step_limit ();
+	}
+      else
+	info(11) = 0;
 
       abs_tol = absolute_tolerance ();
       rel_tol = relative_tolerance ();
@@ -278,65 +300,22 @@
       else
 	{
 	  (*current_liboctave_error_handler)
-	    ("dassl: inconsistent sizes for tolerance arrays");
+	    ("dasrt: inconsistent sizes for tolerance arrays");
 
 	  integration_error = true;
 	  return;
 	}
 
-      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;
-
-      pinfo = info.fortran_vec ();
-      piwork = iwork.fortran_vec ();
       pabs_tol = abs_tol.fortran_vec ();
       prel_tol = rel_tol.fortran_vec ();
-      prwork = rwork.fortran_vec ();
-      pjroot = jroot.fortran_vec ();
 
-      info(5) = 0;
-      info(8) = 0;
-      initialized = true;
+      DASRT_options::reset = false;
     }
 
-  if (restart)
-    {
-      info(0) = 0;
+  static double *dummy = 0;
+  static int *idummy = 0;
 
-      if (stop_time_set)
-	{
-	  info(3) = 1;
-	  rwork(0) = stop_time;
-	}
-      else
-	info(3) = 0;
-    }
-
-  double *dummy = 0;
-  int *idummy = 0;
-
-  F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, px, pxdot, tout, pinfo,
+  F77_XFCN (ddasrt, DASRT, (ddasrt_f, nn, t, px, pxdot, tout, pinfo,
 			    prel_tol, pabs_tol, istate, prwork, lrw,
 			    piwork, liw, dummy, idummy, ddasrt_j,
 			    ddasrt_g, ng, pjroot));
@@ -344,7 +323,7 @@
   if (f77_exception_encountered)
     {
       integration_error = true;
-      (*current_liboctave_error_handler) ("unrecoverable error in dassl");
+      (*current_liboctave_error_handler) ("unrecoverable error in dasrt");
     }
   else
     {
@@ -409,6 +388,7 @@
   ColumnVector t_out = tout;
 
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
@@ -467,6 +447,7 @@
   ColumnVector t_outs = tout;
 
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
--- a/liboctave/DASRT.h
+++ b/liboctave/DASRT.h
@@ -74,12 +74,14 @@
 {
 public:
 
-  DASRT (void);
+  DASRT (void) : DAERT (), DASRT_options (), initialized (false) { }
 
-  DASRT (const ColumnVector& state, double time, DAERTFunc& f);
+  DASRT (const ColumnVector& state, double time, DAERTFunc& f)
+    : DAERT (state, time, f), DASRT_options (), initialized (false) { }
 
   DASRT (const ColumnVector& state, const ColumnVector& deriv,
-	 double time, DAERTFunc& f);
+	 double time, DAERTFunc& f)
+    : DAERT (state, deriv, time, f), DASRT_options (), initialized (false) { }
 
   ~DASRT (void) { }
 
@@ -94,12 +96,9 @@
 
   bool initialized;
 
-  bool sanity_checked;
-
   int liw;  
   int lrw;
 
-  int n;
   int ng;
 
   Array<int> info;
--- a/liboctave/DASSL-opts.in
+++ b/liboctave/DASSL-opts.in
@@ -15,10 +15,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
@@ -35,10 +36,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
--- a/liboctave/DASSL.cc
+++ b/liboctave/DASSL.cc
@@ -56,53 +56,7 @@
 
 static int nn;
 
-DASSL::DASSL (void) : DAE ()
-{
-  liw = 0;
-  lrw = 0;
-
-  sanity_checked = false;
-
-  info.resize (15);
-
-  for (int i = 0; i < 15; i++)
-    info.elem (i) = 0;
-}
-
-DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f)
-  : DAE (state, time, f)
-{
-  n = size ();
-
-  liw = 20 + n;
-  lrw = 40 + 9*n + n*n;
-
-  sanity_checked = false;
-
-  info.resize (15, 0);
-}
-
-DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv,
-	      double time, DAEFunc& f)
-  : DAE (state, deriv, time, f)
-{
-  n = size ();
-
-  DAEFunc::set_function (f.function ());
-  DAEFunc::set_jacobian_function (f.jacobian_function ());
-
-  liw = 20 + n;
-  lrw = 40 + 9*n + n*n;
-
-  sanity_checked = false;
-
-  info.resize (15);
-
-  for (int i = 0; i < 15; i++)
-    info.elem (i) = 0;
-}
-
-int
+static int
 ddassl_f (const double& time, const double *state, const double *deriv,
 	  double *delta, int& ires, double *, int *)
 {
@@ -134,7 +88,7 @@
   return 0;
 }
 
-int
+static int
 ddassl_j (const double& time, const double *state, const double *deriv,
 	  double *pd, const double& cj, double *, int *)
 {
@@ -161,140 +115,158 @@
 ColumnVector
 DASSL::do_integrate (double tout)
 {
-  // XXX FIXME XXX -- should handle all this option stuff just once
-  // for each new problem.
-
   ColumnVector retval;
 
-  if (restart)
+  if (! initialized || restart || DAEFunc::reset|| DASSL_options::reset)
     {
-      restart = false;
-      info.elem (0) = 0;
-    }
+      integration_error = false;
+
+      initialized = true;
+
+      info.resize (15);
+
+      for (int i = 0; i < 15; i++)
+	info(i) = 0;
 
-  if (iwork.length () != liw)
-    iwork.resize (liw);
+      pinfo = info.fortran_vec ();
+
+      int n = size ();
+
+      liw = 20 + n;
+      lrw = 40 + 9*n + n*n;
 
-  if (rwork.length () != lrw)
-    rwork.resize (lrw);
+      nn = n;
 
-  integration_error = false;
+      iwork.resize (liw);
+      rwork.resize (lrw);
 
-  user_fun = DAEFunc::fun;
-  user_jac = DAEFunc::jac;
+      info(0) = 0;
 
-  if (user_jac)
-    info.elem (4) = 1;
-  else
-    info.elem (4) = 0;
+      if (stop_time_set)
+	{
+	  rwork(0) = stop_time;
+	  info(3) = 1;
+	}
+      else
+	info(3) = 0;
 
-  double *px = x.fortran_vec ();
-  double *pxdot = xdot.fortran_vec ();
+      px = x.fortran_vec ();
+      pxdot = xdot.fortran_vec ();
+
+      piwork = iwork.fortran_vec ();
+      prwork = rwork.fortran_vec ();
+
+      restart = false;
+
+      // DAEFunc
 
-  nn = n;
+      user_fun = DAEFunc::function ();
+      user_jac = DAEFunc::jacobian_function ();
+
+      if (user_fun)
+	{
+	  int ires = 0;
+
+	  ColumnVector res = (*user_fun) (x, xdot, t, ires);
 
-  if (! sanity_checked)
-    {
-      int ires = 0;
+	  if (res.length () != x.length ())
+	    {
+	      (*current_liboctave_error_handler)
+		("dassl: inconsistent sizes for state and residual vectors");
 
-      ColumnVector res = (*user_fun) (x, xdot, t, ires);
-
-      if (res.length () != x.length ())
+	      integration_error = true;
+	      return retval;
+	    }
+	}
+      else
 	{
 	  (*current_liboctave_error_handler)
-	    ("dassl: inconsistent sizes for state and residual vectors");
+	    ("dassl: no user supplied RHS subroutine!");
 
 	  integration_error = true;
 	  return retval;
 	}
 
-      sanity_checked = true;
-    }
+      info(4) = user_jac ? 1 : 0;
   
-  if (stop_time_set)
-    {
-      rwork.elem (0) = stop_time;
-      info.elem (3) = 1;
-    }
-  else
-    info.elem (3) = 0;
+      DAEFunc::reset = false;
 
-  Array<double> abs_tol = absolute_tolerance ();
-  Array<double> rel_tol = relative_tolerance ();
+      // DASSL_options
 
-  int abs_tol_len = abs_tol.length ();
-  int rel_tol_len = rel_tol.length ();
+      double hmax = maximum_step_size ();
+      if (hmax >= 0.0)
+	{
+	  rwork(1) = hmax;
+	  info(6) = 1;
+	}
+      else
+	info(6) = 0;
 
-  if (abs_tol_len == 1 && rel_tol_len == 1)
-    {
-      info.elem (1) = 0;
-    }
-  else if (abs_tol_len == n && rel_tol_len == n)
-    {
-      info.elem (1) = 1;
-    }
-  else
-    {
-      (*current_liboctave_error_handler)
-	("dassl: inconsistent sizes for tolerance arrays");
+      double h0 = initial_step_size ();
+      if (h0 >= 0.0)
+	{
+	  rwork(2) = h0;
+	  info(7) = 1;
+	}
+      else
+	info(7) = 0;
 
-      integration_error = true;
-      return retval;
-    }
-
-  double hmax = maximum_step_size ();
-  if (hmax >= 0.0)
-    {
-      rwork.elem (1) = hmax;
-      info.elem (6) = 1;
-    }
-  else
-    info.elem (6) = 0;
+      int maxord = maximum_order ();
+      if (maxord >= 0)
+	{
+	  if (maxord > 0 && maxord < 6)
+	    {
+	      info(8) = 1;
+	      iwork(2) = maxord;
+	    }
+	  else
+	    {
+	      (*current_liboctave_error_handler)
+		("dassl: invalid value for maximum order");
+	      integration_error = true;
+	      return retval;
+	    }
+	}
 
-  double h0 = initial_step_size ();
-  if (h0 >= 0.0)
-    {
-      rwork.elem (2) = h0;
-      info.elem (7) = 1;
-    }
-  else
-    info.elem (7) = 0;
+      int enc = enforce_nonnegativity_constraints ();
+      info(9) = enc ? 1 : 0;
+
+      int ccic = compute_consistent_initial_condition ();
+      info(10) = ccic ? 1 : 0;
+
+      abs_tol = absolute_tolerance ();
+      rel_tol = relative_tolerance ();
 
-  int maxord = maximum_order ();
-  if (maxord >= 0)
-    {
-      if (maxord > 0 && maxord < 6)
+      int abs_tol_len = abs_tol.length ();
+      int rel_tol_len = rel_tol.length ();
+
+      if (abs_tol_len == 1 && rel_tol_len == 1)
 	{
-	  info(8) = 1;
-	  iwork(2) = maxord;
+	  info(1) = 0;
+	}
+      else if (abs_tol_len == n && rel_tol_len == n)
+	{
+	  info(1) = 1;
 	}
       else
 	{
 	  (*current_liboctave_error_handler)
-	    ("dassl: invalid value for maximum order");
+	    ("dassl: inconsistent sizes for tolerance arrays");
+
 	  integration_error = true;
 	  return retval;
 	}
+
+      pabs_tol = abs_tol.fortran_vec ();
+      prel_tol = rel_tol.fortran_vec ();
+
+      DASSL_options::reset = false;
     }
 
-  int enc = enforce_nonnegativity_constraints ();
-  info (9) = enc ? 1 : 0;
-
-  int ccic = compute_consistent_initial_condition ();
-  info(10) = ccic ? 1 : 0;
-
-  double *dummy = 0;
-  int *idummy = 0;
+  static double *dummy = 0;
+  static int *idummy = 0;
 
-  int *pinfo = info.fortran_vec ();
-  int *piwork = iwork.fortran_vec ();
-  double *prwork = rwork.fortran_vec ();
-  double *pabs_tol = abs_tol.fortran_vec ();
-  double *prel_tol = rel_tol.fortran_vec ();
-
-// again:
-
-  F77_XFCN (ddassl, DDASSL, (ddassl_f, n, t, px, pxdot, tout, pinfo,
+  F77_XFCN (ddassl, DDASSL, (ddassl_f, nn, t, px, pxdot, tout, pinfo,
 			     prel_tol, pabs_tol, istate, prwork, lrw,
 			     piwork, liw, dummy, idummy, ddassl_j));
 
@@ -366,7 +338,9 @@
 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out)
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
@@ -409,7 +383,9 @@
 		  const ColumnVector& tcrit) 
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
--- a/liboctave/DASSL.h
+++ b/liboctave/DASSL.h
@@ -37,12 +37,14 @@
 {
 public:
 
-  DASSL (void);
+  DASSL (void) : DAE (), DASSL_options (), initialized (false) { }
 
-  DASSL (const ColumnVector& state, double time, DAEFunc& f);
+  DASSL (const ColumnVector& state, double time, DAEFunc& f)
+    : DAE (state, time, f), DASSL_options (), initialized (false) { }
 
-  DASSL (const ColumnVector& state, const ColumnVector& xdot,
-	 double time, DAEFunc& f);
+  DASSL (const ColumnVector& state, const ColumnVector& deriv,
+	 double time, DAEFunc& f)
+    : DAE (state, deriv, time, f), DASSL_options (), initialized (false) { }
 
   ~DASSL (void) { }
 
@@ -61,20 +63,26 @@
 
 private:
 
-  int n;
+  bool initialized;
+
   int liw;  
   int lrw;
-  bool sanity_checked;
+
   Array<int> info;
   Array<int> iwork;
+
   Array<double> rwork;
 
-  friend int ddassl_j (double *time, double *state, double *deriv,
-		       double *pd, double *cj, double *rpar, int *ipar);
+  Array<double> abs_tol;
+  Array<double> rel_tol;
 
-  friend int ddassl_f (double *time, double *state, double *deriv,
-		       double *delta, int *ires, double *rpar, int *ipar);
-
+  double *px;
+  double *pxdot;
+  double *pabs_tol;
+  double *prel_tol;
+  int *pinfo;
+  int *piwork;
+  double *prwork;
 };
 
 #endif
--- a/liboctave/LSODE-opts.in
+++ b/liboctave/LSODE-opts.in
@@ -15,10 +15,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
--- a/liboctave/LSODE.cc
+++ b/liboctave/LSODE.cc
@@ -54,28 +54,7 @@
 static ODEFunc::ODEJacFunc user_jac;
 static ColumnVector *tmp_x;
 
-LSODE::LSODE (void) : ODE (), LSODE_options ()
-{
-  n = size ();
-
-  itask = 1;
-  iopt = 0;
-
-  sanity_checked = false;
-}
-
-LSODE::LSODE (const ColumnVector& state, double time, const ODEFunc& f)
-  : ODE (state, time, f), LSODE_options ()
-{
-  n = size ();
-
-  itask = 1;
-  iopt = 0;
-
-  sanity_checked = false;
-}
-
-int
+static int
 lsode_f (const int& neq, const double& time, double *,
 	 double *deriv, int& ierr) 
 {
@@ -98,7 +77,7 @@
   return 0;
 }
 
-int
+static int
 lsode_j (const int& neq, const double& time, double *,
 	 const int&, const int&, double *pd, const int& nrowpd)
 {
@@ -122,60 +101,77 @@
 {
   ColumnVector retval;
 
-  if (restart)
+  static int nn = 0;
+
+  if (! initialized || restart || ODEFunc::reset || LSODE_options::reset)
     {
-      restart = false;
-      istate = 1;
-    }
+      integration_error = false;
 
-  if (integration_method () == "stiff")
-    {
-      if (jac)
-	method_flag = 21;
-      else
-	method_flag = 22;
+      initialized = true;
+
+      istate = 1;
+
+      int n = size ();
+
+      nn = n;
 
-      liw = 20 + n;
-      lrw = 22 + n * (9 + n);
-    }
-  else
-    {
-      method_flag = 10;
+      if (integration_method () == "stiff")
+	{
+	  if (jac)
+	    method_flag = 21;
+	  else
+	    method_flag = 22;
 
-      liw = 20;
-      lrw = 22 + 16 * n;
-    }
+	  liw = 20 + n;
+	  lrw = 22 + n * (9 + n);
+	}
+      else
+	{
+	  method_flag = 10;
 
-  if (iwork.length () != liw)
-    {
+	  liw = 20;
+	  lrw = 22 + 16 * n;
+	}
+
       iwork.resize (liw);
 
       for (int i = 4; i < 9; i++)
-	iwork.elem (i) = 0;
-    }
+	iwork(i) = 0;
 
-  if (rwork.length () != lrw)
-    {
       rwork.resize (lrw);
 
       for (int i = 4; i < 9; i++)
-	rwork.elem (i) = 0;
-    }
+	rwork(i) = 0;
 
-  integration_error = false;
+      if (stop_time_set)
+	{
+	  itask = 4;
+	  rwork(0) = stop_time;
+	  iopt = 1;
+	}
+      else
+	{
+	  itask = 1;
+	}
 
-  double *xp = x.fortran_vec ();
+      px = x.fortran_vec ();
 
-  // 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.
+      piwork = iwork.fortran_vec ();
+      prwork = rwork.fortran_vec ();
+
+      restart = false;
+
+      // ODEFunc
 
-  tmp_x = &x;
-  user_fun = function ();
-  user_jac = jacobian_function ();
+      // 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)
-    {
+      tmp_x = &x;
+
+      user_fun = function ();
+      user_jac = jacobian_function ();
+
       ColumnVector xdot = (*user_fun) (x, t);
 
       if (x.length () != xdot.length ())
@@ -187,70 +183,62 @@
 	  return retval;
 	}
 
-      sanity_checked = true;
-    }
+      ODEFunc::reset = false;
+
+      // LSODE_options
+
+      rel_tol = relative_tolerance ();
+      abs_tol = absolute_tolerance ();
+
+      int abs_tol_len = abs_tol.length ();
+
+      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 = true;
+	  return retval;
+	}
 
-  if (stop_time_set)
-    {
-      itask = 4;
-      rwork.elem (0) = stop_time;
-      iopt = 1;
-    }
-  else
-    {
-      itask = 1;
+      double iss = initial_step_size ();
+      if (iss >= 0.0)
+	{
+	  rwork(4) = iss;
+	  iopt = 1;
+	}
+
+      double maxss = maximum_step_size ();
+      if (maxss >= 0.0)
+	{
+	  rwork(5) = maxss;
+	  iopt = 1;
+	}
+
+      double minss = minimum_step_size ();
+      if (minss >= 0.0)
+	{
+	  rwork(6) = minss;
+	  iopt = 1;
+	}
+
+      int sl = step_limit ();
+      if (sl > 0)
+	{
+	  iwork(5) = sl;
+	  iopt = 1;
+	}
+
+      pabs_tol = abs_tol.fortran_vec ();
+
+      LSODE_options::reset = false;
     }
 
-  double rel_tol = relative_tolerance ();
-
-  const Array<double> abs_tol = absolute_tolerance ();
-
-  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 = true;
-      return retval;
-    }
-
-  if (initial_step_size () >= 0.0)
-    {
-      rwork.elem (4) = initial_step_size ();
-      iopt = 1;
-    }
-
-  if (maximum_step_size () >= 0.0)
-    {
-      rwork.elem (5) = maximum_step_size ();
-      iopt = 1;
-    }
-
-  if (minimum_step_size () >= 0.0)
-    {
-      rwork.elem (6) = minimum_step_size ();
-      iopt = 1;
-    }
-
-  if (step_limit () > 0)
-    {
-      iwork.elem (5) = step_limit ();
-      iopt = 1;
-    }
-
-  const double *pabs_tol = abs_tol.fortran_vec ();
-  int *piwork = iwork.fortran_vec ();
-  double *prwork = rwork.fortran_vec ();
-
-  F77_XFCN (lsode, LSODE, (lsode_f, n, xp, t, tout, itol, rel_tol,
+  F77_XFCN (lsode, LSODE, (lsode_f, nn, px, t, tout, itol, rel_tol,
 			   pabs_tol, itask, istate, iopt, prwork, lrw,
 			   piwork, liw, lsode_j, method_flag));
 
@@ -365,7 +353,9 @@
 LSODE::do_integrate (const ColumnVector& tout)
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
@@ -393,7 +383,9 @@
 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
--- a/liboctave/LSODE.h
+++ b/liboctave/LSODE.h
@@ -37,11 +37,10 @@
 {
 public:
 
-  LSODE (void);
+  LSODE (void) : ODE (), LSODE_options (), initialized (false) { }
 
-  LSODE (int n);
-  
-  LSODE (const ColumnVector& state, double time, const ODEFunc& f);
+  LSODE (const ColumnVector& state, double time, const ODEFunc& f)
+    : ODE (state, time, f), LSODE_options (), initialized (false) { }
 
   ~LSODE (void) { }
 
@@ -55,20 +54,27 @@
 
 private:
 
-  int n;
+  bool initialized;
+
   int method_flag;
-  Array<int> iwork;
-  Array<double> rwork;
   int itask;
   int iopt;
+  int itol;
+
   int liw;
   int lrw;
-  bool sanity_checked;
+
+  Array<int> iwork;
+  Array<double> rwork;
+
+  double rel_tol;
 
-  friend int lsode_f (int *neq, double *t, double *y, double *ydot);
+  Array<double> abs_tol;
 
-  friend int lsode_j (int *neq, double *t, double *y, int *ml, int *mu,
-		      double *pd, int *nrowpd);
+  double *px;
+  double *pabs_tol;
+  int *piwork;
+  double *prwork;
 };
 
 #endif
--- a/liboctave/ODEFunc.h
+++ b/liboctave/ODEFunc.h
@@ -35,16 +35,16 @@
   typedef Matrix (*ODEJacFunc) (const ColumnVector&, double);
 
   ODEFunc (void)
-    : fun (0), jac (0) { }
+    : fun (0), jac (0), reset (true) { }
 
   ODEFunc (ODERHSFunc f)
-    : fun (f), jac (0) { }
+    : fun (f), jac (0), reset (true) { }
 
   ODEFunc (ODERHSFunc f, ODEJacFunc j)
-    : fun (f), jac (j) { }
+    : fun (f), jac (j), reset (true) { }
 
   ODEFunc (const ODEFunc& a)
-    : fun (a.fun), jac (a.jac) { }
+    : fun (a.fun), jac (a.jac), reset (true) { }
 
   ODEFunc& operator = (const ODEFunc& a)
     {
@@ -52,6 +52,7 @@
 	{
 	  fun = a.fun;
 	  jac = a.jac;
+	  reset = a.reset;
 	}
       return *this;
     }
@@ -63,6 +64,7 @@
   ODEFunc& set_function (ODERHSFunc f)
     {
       fun = f;
+      reset = true;
       return *this;
     }
 
@@ -71,6 +73,7 @@
   ODEFunc& set_jacobian_function (ODEJacFunc j)
     {
       jac = j;
+      reset = true;
       return *this;
     }
 
@@ -78,6 +81,13 @@
 
   ODERHSFunc fun;
   ODEJacFunc jac;
+
+  // This variable is TRUE when this object is constructed, and also
+  // after any internal data has changed.  Derived classes may use
+  // this information (and change it) to know when to (re)initialize
+  // their own internal data related to this object.
+
+  bool reset;
 };
 
 #endif
--- a/liboctave/ODESSA-opts.in
+++ b/liboctave/ODESSA-opts.in
@@ -15,10 +15,11 @@
       {
         $OPTVAR.resize (1);
         $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON);
+        reset = true;
       }
 
     void set_$OPT (const $TYPE& val)
-      { $OPTVAR = val; }
+      { $OPTVAR = val; reset = true; }
   END_SET_CODE
 END_OPTION
 
--- a/liboctave/ODESSA.h
+++ b/liboctave/ODESSA.h
@@ -109,14 +109,13 @@
   int n;
   int npar;
 
-  // Hey, check out this crap: ZZZZ
+  // XXX FIXME XXX -- ???
   Array<double> par;
 
   Matrix sx0;
 
   Matrix y;
 
-
   double *py;
   double *ppar;
   int *piwork;
--- a/liboctave/base-de.h
+++ b/liboctave/base-de.h
@@ -82,9 +82,14 @@
     {
       stop_time_set = true;
       stop_time = t;
+      force_restart ();
     }
 
-  void clear_stop_time (void) { stop_time_set = false; }
+  void clear_stop_time (void)
+    {
+      stop_time_set = false;
+      force_restart ();
+    }
 
   virtual void force_restart (void) { restart = true; }
 
--- a/mk-opts.pl
+++ b/mk-opts.pl
@@ -374,7 +374,8 @@
         }
     }
 
-  print "    }\n";
+  print "      reset = true;
+    }\n";
 
   print "\n  void copy (const ${class_name}& opt)\n    {\n";
 
@@ -383,7 +384,8 @@
       print "      $optvar[$i] = opt.$optvar[$i];\n";
     }
 
-  print "    }\n";
+  print "      reset = opt.reset;
+    }\n";
 
   print "\n  void set_default_options (void) { init (); }\n";
 
@@ -393,7 +395,7 @@
         {
           &emit_set_decl ($i);
 
-          print "\n    { $optvar[$i] = $set_expr[$i]; }\n";
+          print "\n    { $optvar[$i] = $set_expr[$i]; reset = true; }\n";
         }
       elsif ($set_body[$i])
         {
@@ -403,7 +405,7 @@
           chop ($s);
           $s =~ s/^/  /g;
           $s =~ s/\n/\n  /g;
-          print "\n    {\n$s\n    }\n";
+          print "\n    {\n$s\n      reset = true;\n    }\n";
         }
       elsif ($set_code[$i])
         {
@@ -427,7 +429,7 @@
       print "  $type[$i] $optvar[$i];\n";
     }
 
-  print "};\n\n#endif\n";
+  print "\nprotected:\n\n  bool reset;\n};\n\n#endif\n";
 }
 
 sub emit_set_decl