--- a/libinterp/dldfcn/__ode15__.cc	Fri Dec 24 18:51:53 2021 +0300
+++ b/libinterp/dldfcn/__ode15__.cc	Sat Dec 25 12:36:11 2021 +0100
@@ -194,6 +194,9 @@
       IDAFree (&m_mem);
       SUNLinSolFree (m_sunLinearSolver);
       SUNMatDestroy (m_sunJacMatrix);
+#  if defined (HAVE_SUNDIALS_SUNCONTEXT)
+      SUNContext_Free (&m_sunContext);
+#  endif
     }
 
     IDA&
@@ -253,7 +256,11 @@
 
     static ColumnVector NVecToCol (N_Vector& v, octave_f77_int_type n);
 
+#  if defined (HAVE_SUNDIALS_SUNCONTEXT)
+    N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
+#  else
     static N_Vector ColToNVec (const ColumnVector& data, octave_f77_int_type n);
+#  endif
 
     void
     set_up (const ColumnVector& y);
@@ -366,6 +373,9 @@
     DAEJacFuncSparse m_jacspfun;
     DAEJacCellDense m_jacdcell;
     DAEJacCellSparse m_jacspcell;
+#  if defined (HAVE_SUNDIALS_SUNCONTEXT)
+    SUNContext m_sunContext;
+#  endif
     SUNMatrix m_sunJacMatrix;
     SUNLinearSolver m_sunLinearSolver;
   };
@@ -395,6 +405,12 @@
       puntrr[i] = res(i);
   }
 
+#  if defined (HAVE_SUNDIALS_SUNCONTEXT)
+#    define OCTAVE_SUNCONTEXT , m_sunContext
+#  else
+#    define OCTAVE_SUNCONTEXT
+#  endif
+
   void
   IDA::set_up (const ColumnVector& y)
   {
@@ -407,18 +423,21 @@
         // Initially allocate memory for 0 entries. We will reallocate when we
         // get the Jacobian matrix from the user and know the actual number of
         // entries.
-        m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT);
+        m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, 0, CSC_MAT
+                                          OCTAVE_SUNCONTEXT);
 #    else
         octave_f77_int_type max_elems;
         if (math::int_multiply_overflow (m_num, m_num, &max_elems))
           error ("Unable to allocate memory for sparse Jacobian");
 
-        m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, max_elems, CSC_MAT);
+        m_sunJacMatrix = SUNSparseMatrix (m_num, m_num, max_elems, CSC_MAT
+                                          OCTAVE_SUNCONTEXT);
 #    endif
         if (! m_sunJacMatrix)
           error ("Unable to create sparse Jacobian for Sundials");
 
-        m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix);
+        m_sunLinearSolver = SUNLinSol_KLU (yy, m_sunJacMatrix
+                                           OCTAVE_SUNCONTEXT);
         if (! m_sunLinearSolver)
           error ("Unable to create KLU sparse solver");
 
@@ -437,11 +456,12 @@
     else
       {
 
-        m_sunJacMatrix = SUNDenseMatrix (m_num, m_num);
+        m_sunJacMatrix = SUNDenseMatrix (m_num, m_num OCTAVE_SUNCONTEXT);
         if (! m_sunJacMatrix)
           error ("Unable to create dense Jacobian for Sundials");
 
-        m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix);
+        m_sunLinearSolver = SUNLinSol_Dense (yy, m_sunJacMatrix
+                                             OCTAVE_SUNCONTEXT);
         if (! m_sunLinearSolver)
           error ("Unable to create dense linear solver");
 
@@ -540,7 +560,7 @@
   N_Vector
   IDA::ColToNVec (const ColumnVector& data, octave_f77_int_type n)
   {
-    N_Vector v = N_VNew_Serial (n);
+    N_Vector v = N_VNew_Serial (n OCTAVE_SUNCONTEXT);
 
     realtype *punt = nv_data_s (v);
 
@@ -563,7 +583,13 @@
   IDA::initialize (void)
   {
     m_num = to_f77_int (m_y0.numel ());
+#  if defined (HAVE_SUNDIALS_SUNCONTEXT)
+    if (SUNContext_Create (nullptr, &m_sunContext) < 0)
+      error ("__ode15__: unable to create context for SUNDIALS");
+    m_mem = IDACreate (m_sunContext);
+#  else
     m_mem = IDACreate ();
+#  endif
 
     N_Vector yy = ColToNVec (m_y0, m_num);
 
@@ -726,7 +752,7 @@
         if (status == 0)
           {
             // Interpolate in tend
-            N_Vector dky = N_VNew_Serial (m_num);
+            N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
 
             if (IDAGetDky (m_mem, tend, 0, dky) != 0)
               error ("IDA failed to interpolate y");
@@ -890,9 +916,9 @@
     realtype h = 0, tcur = 0;
     bool status = 0;
 
-    N_Vector dky = N_VNew_Serial (m_num);
+    N_Vector dky = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
 
-    N_Vector dkyp = N_VNew_Serial (m_num);
+    N_Vector dkyp = N_VNew_Serial (m_num OCTAVE_SUNCONTEXT);
 
     ColumnVector yout (m_num);
     ColumnVector ypout (m_num);

--- a/m4/acinclude.m4	Fri Dec 24 18:51:53 2021 +0300
+++ b/m4/acinclude.m4	Sat Dec 25 12:36:11 2021 +0100
@@ -2378,13 +2378,9 @@
     octave_have_sundials_compatible_api=no
   fi
   AC_MSG_RESULT([$octave_have_sundials_compatible_api])
-  dnl Octave doesn't yet support the SUNContext API introduced in SUNDIALS 6.0.
-  dnl For now, check for that API and de-activate features if it is found.
-  dnl FIXME: Properly support that API.
-  AC_MSG_CHECKING([whether SUNDIALS API uses SUNContext object])
-  AC_MSG_RESULT([$ac_cv_func_SUNContext_Create])
   if test "x$ac_cv_func_SUNContext_Create" = xyes; then
-    octave_have_sundials_compatible_api=no
+    AC_DEFINE(HAVE_SUNDIALS_SUNCONTEXT, 1,
+      [Define to 1 if SUNDIALS' API is using a SUNContext object.])
   fi
   if test $octave_have_sundials_compatible_api = no; then
     warn_sundials_disabled="SUNDIALS libraries do not provide an API that is compatible with Octave.  The solvers ode15i and ode15s will be disabled."
@@ -2452,11 +2448,12 @@
      #  include <ufsparse/klu.h>
      #endif
     ])
+  ## Check for current KLU function name first.
   OCTAVE_CHECK_LIB(sundials_sunlinsolklu, SUNLINSOL_KLU, [],
-    [], [SUNKLU], [],
+    [], [SUNLinSol_KLU], [],
     [don't use SUNDIALS SUNLINSOL_KLU library, disable ode15i and ode15s sparse Jacobian],
-    [AC_CHECK_FUNCS([SUNLinSol_KLU SUNKLU])
-     AC_CACHE_CHECK([whether compiling a program that calls SUNKLU works],
+    [AC_CHECK_FUNCS([SUNLinSol_KLU])
+     AC_CACHE_CHECK([whether compiling a program that calls SUNLinSol_KLU works],
       [octave_cv_sundials_sunlinsol_klu],
       [AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
          #if defined (HAVE_IDA_IDA_H)
@@ -2478,11 +2475,53 @@
          #include <sunlinsol/sunlinsol_klu.h>
          #endif
          ]], [[
-         SUNKLU (0, 0);
+         #if defined (HAVE_SUNCONTEXT_CREATE)
+           SUNContext *sunContext;
+           if (SUNContext_Create (NULL, sunContext) < 0)
+             1/0;  // provoke an error
+           SUNLinSol_KLU (0, 0, *sunContext);
+           SUNContext_Free (sunContext);
+         #else
+           SUNLinSol_KLU (0, 0);
+         #endif
       ]])],
       octave_cv_sundials_sunlinsol_klu=yes,
       octave_cv_sundials_sunlinsol_klu=no)
     ])])
+  if test "x$octave_cv_sundials_sunlinsol_klu" = xno; then
+    ## Check for deprecated KLU function name second.
+    OCTAVE_CHECK_LIB(sundials_sunlinsolklu, SUNLINSOL_KLU, [],
+      [], [SUNKLU], [],
+      [don't use SUNDIALS SUNLINSOL_KLU library, disable ode15i and ode15s sparse Jacobian],
+      [AC_CHECK_FUNCS([SUNKLU])
+       AC_CACHE_CHECK([whether compiling a program that calls SUNLinSol_KLU works],
+        [octave_cv_sundials_sunlinsol_klu],
+        [AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
+           #if defined (HAVE_IDA_IDA_H)
+           #include <ida/ida.h>
+           #endif
+           #if defined (HAVE_KLU_H)
+           #include <klu.h>
+           #endif
+           #if defined (HAVE_KLU_KLU_H)
+           #include <klu/klu.h>
+           #endif
+           #if defined (HAVE_SUITESPARSE_KLU_H)
+           #include <suitesparse/klu.h>
+           #endif
+           #if defined (HAVE_UFPARSE_KLU_H)
+           #include <ufsparse/klu.h>
+           #endif
+           #if defined (HAVE_SUNLINSOL_SUNLINSOL_KLU_H)
+           #include <sunlinsol/sunlinsol_klu.h>
+           #endif
+           ]], [[
+           SUNKLU (0, 0);
+        ]])],
+        octave_cv_sundials_sunlinsol_klu=yes,
+        octave_cv_sundials_sunlinsol_klu=no)
+      ])])
+  fi
   if test "x$ac_cv_header_sunlinsol_sunlinsol_klu_h" = xyes \
      && test "x$octave_cv_sundials_sunlinsol_klu" = xyes; then
     AC_DEFINE(HAVE_SUNDIALS_SUNLINSOL_KLU, 1,