--- a/build-aux/mk-opts.pl	Thu Sep 29 23:09:05 2022 -0400
+++ b/build-aux/mk-opts.pl	Fri Nov 12 08:53:05 2010 +0100
@@ -427,7 +427,7 @@
 
   for ($i = 0; $i < $OPT_NUM; $i++)
     {
-      if ($INIT_VALUE[$i])
+      if (defined $INIT_VALUE[$i])
         {
           print "      $OPTVAR[$i] = $INIT_VALUE[$i];\n";
         }

--- a/liboctave/numeric/LSODE-opts.in	Thu Sep 29 23:09:05 2022 -0400
+++ b/liboctave/numeric/LSODE-opts.in	Fri Nov 12 08:53:05 2010 +0100
@@ -164,3 +164,67 @@
   INIT_VALUE = "100000"
   SET_EXPR = "val"
 END_OPTION
+
+OPTION
+  NAME = "jacobian type"
+  DOC_ITEM
+A string specifying the type of Jacobian used with the stiff backward
+differentiation formula (BDF) integration method. Valid values are
+
+@table @asis
+@item @qcode{"full"}
+The default. All partial derivatives are approximated or used from the
+user-supplied Jacobian function.
+
+@item @qcode{"banded"}
+Only the diagonal and the number of lower and upper subdiagonals specified by
+the options @qcode{"lower jacobian subdiagonals"} and @qcode{"upper jacobian
+subdiagonals"}, respectively, are approximated or used from the user-supplied
+Jacobian function. A user-supplied Jacobian function may set all other
+partial derivatives to arbitrary values.
+
+@item @qcode{"diagonal"}
+If a Jacobian function is supplied by the user, this setting has no effect.
+A Jacobian approximated by @code{lsode} is restricted to the diagonal, where
+each partial derivative is computed by applying a finite change to all
+elements of the state together; if the real Jacobian is indeed always diagonal,
+this has the same effect as applying the finite change only to the respective
+element of the state, but is more efficient.
+@end table
+
+  END_DOC_ITEM
+  TYPE = "std::string"
+  SET_ARG_TYPE = "const $TYPE&"
+  INIT_VALUE = ""full""
+  SET_BODY
+    if (val == "full" || val == "banded" || val == "diagonal")
+      $OPTVAR = val;
+    else
+      (*current_liboctave_error_handler)
+        ("lsode_options: jacobian type must be \"full\", \"banded\", or \"diagonal\"");
+  END_SET_BODY
+END_OPTION
+
+OPTION
+  NAME = "lower jacobian subdiagonals"
+  DOC_ITEM
+Number of lower subdiagonals used if option @qcode{"jacobian type"} is set to
+@qcode{"banded"}. The default is zero.
+
+  END_DOC_ITEM
+  TYPE = "octave_idx_type"
+  INIT_VALUE = "0"
+  SET_EXPR = "val"
+END_OPTION
+
+OPTION
+  NAME = "upper jacobian subdiagonals"
+  DOC_ITEM
+Number of upper subdiagonals used if option @qcode{"jacobian type"} is set to
+@qcode{"banded"}. The default is zero.
+
+  END_DOC_ITEM
+  TYPE = "octave_idx_type"
+  INIT_VALUE = "0"
+  SET_EXPR = "val"
+END_OPTION

--- a/liboctave/numeric/LSODE.cc	Thu Sep 29 23:09:05 2022 -0400
+++ b/liboctave/numeric/LSODE.cc	Fri Nov 12 08:53:05 2010 +0100
@@ -55,6 +55,7 @@
 static ODEFunc::ODERHSFunc user_fcn;
 static ODEFunc::ODEJacFunc user_jac;
 static ColumnVector *tmp_x;
+static bool user_jac_ignore_ml_mu;
 
 static F77_INT
 lsode_f (const F77_INT& neq, const double& time, double *, double *deriv,
@@ -80,8 +81,9 @@
 }
 
 static F77_INT
-lsode_j (const F77_INT& neq, const double& time, double *, const F77_INT&,
-         const F77_INT&, double *pd, const F77_INT& nrowpd)
+lsode_j (const F77_INT& neq, const double& time, double *, 
+         const F77_INT& ml, const F77_INT& mu,
+         double *pd, const F77_INT& nrowpd)
 {
   Matrix tmp_jac (neq, neq);
 
@@ -91,10 +93,17 @@
 
   tmp_jac = (*user_jac) (*tmp_x, time);
 
-  for (F77_INT j = 0; j < neq; j++)
-    for (F77_INT i = 0; i < neq; i++)
-      pd[nrowpd * j + i] = tmp_jac (i, j);
-
+  if (user_jac_ignore_ml_mu)
+    for (F77_INT j = 0; j < neq; j++)
+      for (F77_INT i = 0; i < neq; i++)
+        pd[nrowpd * j + i] = tmp_jac (i, j);
+  else
+    // upper left ends of subdiagonals in tmp_jac
+    for (F77_INT i = 0, j = mu; i <= ml; j == 0 ? i++ : j--)
+      // walk down the subdiagonal in tmp_jac
+      for (F77_INT k = i, l = j; k < neq && l < neq; k++, l++)
+        pd[nrowpd * l + k + mu - l] = tmp_jac (k, l);    
+      
   return 0;
 }
 
@@ -120,14 +129,43 @@
 
       octave_idx_type max_maxord = 0;
 
+      user_jac_ignore_ml_mu = true;
+
+      m_iwork(0) = lower_jacobian_subdiagonals (); // 'ML' in dlsode.f
+
+      m_iwork(1) = upper_jacobian_subdiagonals (); // 'MU' in dlsode.f
+
       if (integration_method () == "stiff")
         {
           max_maxord = 5;
 
           if (m_jac)
-            m_method_flag = 21;
+            {
+              if (jacobian_type () == "banded")
+                {
+                  m_method_flag = 24;
+                  user_jac_ignore_ml_mu = false;
+                }
+              else
+                m_method_flag = 21;
+            }
           else
-            m_method_flag = 22;
+            {
+              if (jacobian_type () == "full")
+                m_method_flag = 22;
+              else if (jacobian_type () == "banded")
+                m_method_flag = 25;
+              else if (jacobian_type () == "diagonal")
+                m_method_flag = 23;
+              else
+                {
+                  // should be prevented by lsode_options
+                  (*current_liboctave_error_handler)
+                    ("lsode: internal error, wrong jacobian type");
+                  m_integration_error = true;
+                  return retval;
+                }
+            }
 
           m_liw = 20 + n;
           m_lrw = 22 + n * (9 + n);