--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,5 +1,9 @@
 2002-07-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* DAEFunc.h (DAEFunc): Jacobian function now follows format of DASSL.
+	* DASSL.cc (ddassl_j): Make it work.
+	* DASPK.cc (ddaspk_j): Likewise.
+
 	* DAE.cc: Delete.
 
 	* DAERT.h, DAERTFunc.h, DASRT.h, DASRT.cc: New files for DAE

--- a/liboctave/DAEFunc.h
+++ b/liboctave/DAEFunc.h
@@ -31,18 +31,17 @@
 {
 public:
 
-  struct DAEJac
-    {
-      Matrix *dfdxdot;
-      Matrix *dfdx;
-    };
-
   typedef ColumnVector (*DAERHSFunc) (const ColumnVector& x,
 				      const ColumnVector& xdot,
 				      double t, int& ires);
 
-  typedef DAEJac (*DAEJacFunc) (const ColumnVector& x,
-				const ColumnVector& xdot, double t);
+  // This is really the form used by DASSL:
+  //
+  //   PD = DG/DY + CJ * DG/DYPRIME
+
+  typedef Matrix (*DAEJacFunc) (const ColumnVector& x,
+				const ColumnVector& xdot,
+				double t, double cj);
 
   DAEFunc (void)
     : fun (0), jac (0) { }

--- a/liboctave/DASPK.cc
+++ b/liboctave/DASPK.cc
@@ -175,31 +175,27 @@
   return 0;
 }
 
+
 int
-ddaspk_j (const double& time, const double *, const double *, double *pd,
-	  const double& cj, double *, int *)
+ddaspk_j (const double& time, const double *state, const double *deriv,
+	  double *pd, const double& cj, double *, int *)
 {
+  // XXX FIXME XXX -- would be nice to avoid copying the data.
+
   ColumnVector tmp_state (nn);
   ColumnVector tmp_deriv (nn);
 
-  // XXX FIXME XXX
-
-  Matrix tmp_dfdxdot (nn, nn);
-  Matrix tmp_dfdx (nn, nn);
+  for (int i = 0; i < nn; i++)
+    {
+      tmp_deriv.elem (i) = deriv [i];
+      tmp_state.elem (i) = state [i];
+    }
 
-  DAEFunc::DAEJac tmp_jac;
-  tmp_jac.dfdxdot = &tmp_dfdxdot;
-  tmp_jac.dfdx    = &tmp_dfdx;
-
-  tmp_jac = user_jac (tmp_state, tmp_deriv, time);
-
-  // Fix up the matrix of partial derivatives for daspk.
-
-  tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot;
+  Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj);
 
   for (int j = 0; j < nn; j++)
     for (int i = 0; i < nn; i++)
-      pd [nn * j + i] = tmp_dfdx.elem (i, j);
+      pd [nn * j + i] = tmp_pd.elem (i, j);
 
   return 0;
 }

--- a/liboctave/DASRT.cc
+++ b/liboctave/DASRT.cc
@@ -45,10 +45,6 @@
 #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"
@@ -57,23 +53,22 @@
 #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*);
+
+typedef int (*dasrt_jac_ptr) (const double&, const double*, const double*,
+			      double*, const double&, double*, int*);
+
+typedef int (*dasrt_constr_ptr) (const int&, const double&, const double*,
+				 const int&, double*, double*, int*);
+
 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>;
+int F77_FUNC (ddasrt, DASRT) (dasrt_fcn_ptr, const int&, double&,
+			      double*, double*, const double&, int*,
+			      double*, double*, int&, double*,
+			      const int&, int*, const int&, double*,
+			      int*, dasrt_jac_ptr, dasrt_constr_ptr,
+			      const int&, int*);
 
 static DAEFunc::DAERHSFunc user_fsub;
 static DAEFunc::DAEJacFunc user_jsub;
@@ -81,8 +76,8 @@
 static int nn;
 
 static int
-ddasrt_f (const double& t, double *state, double *deriv, double *delta,
-          int& ires, double *rpar, int *ipar)
+ddasrt_f (const double& t, const double *state, const double *deriv,
+	  double *delta, int& ires, double *rpar, int *ipar)
 {
   ColumnVector tmp_state (nn);
   for (int i = 0; i < nn; i++)
@@ -111,43 +106,33 @@
 
 //static efptr e_fun;
 
-static int
-ddasrt_j (const double& t, double *state, double *deriv,
-	  double *pdwork, const double& cj, double *rpar, int *ipar) 
+int
+ddasrt_j (const double& time, const double *state, const double *deriv,
+	  double *pd, const double& cj, double *, int *)
 {
-  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 -- would be nice to avoid copying the data.
 
-  // XXX FIXME XXX
-
-  Matrix tmp_dfdxdot (nn, nn);
-  Matrix tmp_dfdx (nn, nn);
+  ColumnVector tmp_state (nn);
+  ColumnVector tmp_deriv (nn);
 
-  DAEFunc::DAEJac tmp_jac;
-  tmp_jac.dfdxdot = &tmp_dfdxdot;
-  tmp_jac.dfdx    = &tmp_dfdx;
+  for (int i = 0; i < nn; i++)
+    {
+      tmp_deriv.elem (i) = deriv [i];
+      tmp_state.elem (i) = state [i];
+    }
 
-  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;
+  Matrix tmp_pd = user_jsub (tmp_state, tmp_deriv, time, cj);
 
   for (int j = 0; j < nn; j++)
     for (int i = 0; i < nn; i++)
-      pdwork[j*nn+i] = tmp_dfdx.elem (i, j);
+      pd [nn * j + i] = tmp_pd.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) 
+ddasrt_g (const int& neq, const double& t, const double *state,
+	  const int& ng, double *gout, double *rpar, int *ipar) 
 {
   int n = neq;
 

--- a/liboctave/DASSL.cc
+++ b/liboctave/DASSL.cc
@@ -35,19 +35,19 @@
 #include "f77-fcn.h"
 #include "lo-error.h"
 
-typedef int (*dassl_fcn_ptr) (const double&, double*, double*,
+typedef int (*dassl_fcn_ptr) (const double&, const double*, const double*,
 			      double*, int&, double*, int*);
 
-typedef int (*dassl_jac_ptr) (const double&, double*, double*,
+typedef int (*dassl_jac_ptr) (const double&, const double*, const double*,
 			      double*, const double&, double*, int*);
 
 extern "C"
 int F77_FUNC (ddassl, DDASSL) (dassl_fcn_ptr, const int&, double&,
-			      double*, double*, double&, const int*,
-			      const double&, const double&, int&,
-			      double*, const int&, int*, const int&,
-			      const double*, const int*,
-			      dassl_jac_ptr);
+			       double*, double*, double&, const int*,
+			       const double&, const double&, int&,
+			       double*, const int&, int*, const int&,
+			       const double*, const int*,
+			       dassl_jac_ptr);
 
 static DAEFunc::DAERHSFunc user_fun;
 static DAEFunc::DAEJacFunc user_jac;
@@ -132,9 +132,11 @@
 }
 
 int
-ddassl_f (const double& time, double *state, double *deriv,
+ddassl_f (const double& time, const double *state, const double *deriv,
 	  double *delta, int& ires, double *, int *)
 {
+  // XXX FIXME XXX -- would be nice to avoid copying the data.
+
   ColumnVector tmp_deriv (nn);
   ColumnVector tmp_state (nn);
   ColumnVector tmp_delta (nn);
@@ -162,30 +164,25 @@
 }
 
 int
-ddassl_j (const double& time, double *, double *, double *pd, const
-	  double& cj, double *, int *)
+ddassl_j (const double& time, const double *state, const double *deriv,
+	  double *pd, const double& cj, double *, int *)
 {
+  // XXX FIXME XXX -- would be nice to avoid copying the data.
+
   ColumnVector tmp_state (nn);
   ColumnVector tmp_deriv (nn);
 
-  // XXX FIXME XXX
-
-  Matrix tmp_dfdxdot (nn, nn);
-  Matrix tmp_dfdx (nn, nn);
+  for (int i = 0; i < nn; i++)
+    {
+      tmp_deriv.elem (i) = deriv [i];
+      tmp_state.elem (i) = state [i];
+    }
 
-  DAEFunc::DAEJac tmp_jac;
-  tmp_jac.dfdxdot = &tmp_dfdxdot;
-  tmp_jac.dfdx    = &tmp_dfdx;
-
-  tmp_jac = user_jac (tmp_state, tmp_deriv, time);
-
-  // Fix up the matrix of partial derivatives for dassl.
-
-  tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot;
+  Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj);
 
   for (int j = 0; j < nn; j++)
     for (int i = 0; i < nn; i++)
-      pd [nn * j + i] = tmp_dfdx.elem (i, j);
+      pd [nn * j + i] = tmp_pd.elem (i, j);
 
   return 0;
 }

--- a/liboctave/ODESSA.cc
+++ b/liboctave/ODESSA.cc
@@ -67,14 +67,14 @@
 }
 
 typedef int (*odessa_fcn_ptr) (int*, const double&, double*,
-                              double*, double*);
+			       double*, double*);
 
 typedef int (*odessa_jac_ptr) (int*, const double&, double*,
-                              double*, const int&, const int&,
-                              double*, const int&);
+			       double*, const int&, const int&,
+			       double*, const int&);
 
 typedef int (*odessa_dfdp_ptr) (int*, const double&, double*,
-                              double*, double*, const int&);
+				double*, double*, const int&);
 
 
 extern "C"

--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,5 +1,8 @@
 2002-07-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* DLD-FUNCTIONS/dassl.cc (dassl_user_jacobian): New function.
+	(Fdassl): Handle Jacobian function.
+
 	* DLD-FUNCTIONS/dasrt.cc: New file.
 	* Makefile.in (DLD_XSRC): Add it to the list.
 

--- a/src/DLD-FUNCTIONS/dassl.cc
+++ b/src/DLD-FUNCTIONS/dassl.cc
@@ -44,6 +44,9 @@
 // Global pointer for user defined function required by dassl.
 static octave_function *dassl_fcn;
 
+// Global pointer for optional user defined jacobian function.
+static octave_function *dassl_jac;
+
 static DASSL_options dassl_opts;
 
 // Is this a recursive call?
@@ -114,6 +117,70 @@
   return retval;
 }
 
+Matrix
+dassl_user_jacobian (const ColumnVector& x, const ColumnVector& xdot,
+		     double t, double cj)
+{
+  Matrix retval;
+
+  int nstates = x.capacity ();
+
+  assert (nstates == xdot.capacity ());
+
+  octave_value_list args;
+
+  args(3) = cj;
+  args(2) = t;
+
+  if (nstates > 1)
+    {
+      Matrix m1 (nstates, 1);
+      Matrix m2 (nstates, 1);
+      for (int i = 0; i < nstates; i++)
+	{
+	  m1 (i, 0) = x (i);
+	  m2 (i, 0) = xdot (i);
+	}
+      octave_value state (m1);
+      octave_value deriv (m2);
+      args(1) = deriv;
+      args(0) = state;
+    }
+  else
+    {
+      double d1 = x (0);
+      double d2 = xdot (0);
+      octave_value state (d1);
+      octave_value deriv (d2);
+      args(1) = deriv;
+      args(0) = state;
+    }
+
+  if (dassl_jac)
+    {
+      octave_value_list tmp = dassl_jac->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("dassl");
+	  return retval;
+	}
+
+      int tlen = tmp.length ();
+      if (tlen > 0 && tmp(0).is_defined ())
+	{
+	  retval = tmp(0).matrix_value ();
+
+	  if (error_state || retval.length () == 0)
+	    gripe_user_supplied_eval ("dassl");
+	}
+      else
+	gripe_user_supplied_eval ("dassl");
+    }
+
+  return retval;
+}
+
 #define DASSL_ABORT() \
   do \
     { \
@@ -138,7 +205,7 @@
     } \
   while (0)
 
-DEFUN_DLD (dassl, args, ,
+DEFUN_DLD (dassl, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {[@var{x}, @var{xdot}] =} dassl (@var{fcn}, @var{x0}, @var{xdot0}, @var{t}, @var{t_crit})\n\
 Return a matrix of states and their first derivatives with respect to\n\
@@ -191,14 +258,48 @@
 
   int nargin = args.length ();
 
-  if (nargin > 3 && nargin < 6)
+  if (nargin > 3 && nargin < 6 && nargout < 5)
     {
-      dassl_fcn = extract_function
-	(args(0), "dassl", "__dassl_fcn__",
-	 "function res = __dassl_fcn__ (x, xdot, t) res = ",
-	 "; endfunction");
+      dassl_fcn = 0;
+      dassl_jac = 0;
+
+      octave_value f_arg = args(0);
+
+      switch (f_arg.rows ())
+	{
+	case 1:
+	  dassl_fcn = extract_function
+	    (f_arg, "dassl", "__dassl_fcn__",
+	     "function res = __dassl_fcn__ (x, xdot, t) res = ",
+	     "; endfunction");
+	  break;
+
+	case 2:
+	  {
+	    string_vector tmp = f_arg.all_strings ();
 
-      if (! dassl_fcn)
+	    if (! error_state)
+	      {
+		dassl_fcn = extract_function
+		  (tmp(0), "dassl", "__dassl_fcn__",
+		   "function res = __dassl_fcn__ (x, xdot, t) res = ",
+		   "; endfunction");
+
+		if (dassl_fcn)
+		  {
+		    dassl_jac = extract_function
+		      (tmp(1), "dassl", "__dassl_jac__",
+		       "function jac = __dassl_jac__ (x, xdot, t, cj) jac = ",
+		       "; endfunction");
+
+		    if (! dassl_jac)
+		      dassl_fcn = 0;
+		  }
+	      }
+	  }
+	}
+
+      if (error_state || ! dassl_fcn)
 	DASSL_ABORT ();
 
       ColumnVector state = ColumnVector (args(1).vector_value ());
@@ -234,7 +335,11 @@
       double tzero = out_times (0);
 
       DAEFunc func (dassl_user_function);
+      if (dassl_jac)
+	func.set_jacobian_function (dassl_user_jacobian);
+
       DASSL dae (state, deriv, tzero, func);
+
       dae.copy (dassl_opts);
 
       Matrix output;
@@ -247,10 +352,8 @@
 
       if (! error_state)
 	{
-	  retval.resize (2);
-
+	  retval(1) = deriv_output;
 	  retval(0) = output;
-	  retval(1) = deriv_output;
 	}
     }
   else
@@ -394,7 +497,7 @@
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {} dassl_options (@var{opt}, @var{val})\n\
 When called with two arguments, this function allows you set options\n\
-parameters for the function @code{lsode}.  Given one argument,\n\
+parameters for the function @code{dassl}.  Given one argument,\n\
 @code{dassl_options} returns the value of the corresponding option.  If\n\
 no arguments are supplied, the names of all the available options and\n\
 their current values are displayed.\n\

--- a/src/DLD-FUNCTIONS/lsode.cc
+++ b/src/DLD-FUNCTIONS/lsode.cc
@@ -201,20 +201,23 @@
 
   if (nargin > 2 && nargin < 5 && nargout < 4)
     {
+      lsode_fcn = 0;
+      lsode_jac = 0;
+
       octave_value f_arg = args(0);
 
       switch (f_arg.rows ())
 	{
 	case 1:
 	  lsode_fcn = extract_function
-	    (args(0), "lsode", "__lsode_fcn__",
+	    (f_arg, "lsode", "__lsode_fcn__",
 	     "function xdot = __lsode_fcn__ (x, t) xdot = ",
 	     "; endfunction");
 	  break;
 
 	case 2:
 	  {
-	    string_vector tmp = args(0).all_strings ();
+	    string_vector tmp = f_arg.all_strings ();
 
 	    if (! error_state)
 	      {