diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 1586f18106..a728a307f4 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -468,6 +468,8 @@ jobs:
           fmt-ver: 9.1
         - sundials-ver: 6.4.1
           fmt-ver: 10
+        - sundials-ver: 6.6
+          fmt-ver: 10
       fail-fast: false
     steps:
       - uses: actions/checkout@v3
diff --git a/SConstruct b/SConstruct
index 7d9e2d1fdf..3555061271 100644
--- a/SConstruct
+++ b/SConstruct
@@ -1512,7 +1512,7 @@ if env['system_sundials'] == 'y':
     if sundials_ver < parse_version("3.0") or sundials_ver >= parse_version("7.0"):
         logger.error(f"Sundials version {env['sundials_version']!r} is not supported.")
         sys.exit(1)
-    elif sundials_ver > parse_version("6.4.1"):
+    elif sundials_ver > parse_version("6.6.0"):
         logger.warning(f"Sundials version {env['sundials_version']!r} has not been tested.")
 
     logger.info(f"Using system installation of Sundials version {sundials_version!r}.")
diff --git a/include/cantera/kinetics/ImplicitSurfChem.h b/include/cantera/kinetics/ImplicitSurfChem.h
index d4fbd31360..9d8de3945d 100644
--- a/include/cantera/kinetics/ImplicitSurfChem.h
+++ b/include/cantera/kinetics/ImplicitSurfChem.h
@@ -150,7 +150,7 @@ class ImplicitSurfChem : public FuncEval
     // overloaded methods of class FuncEval
 
     //! Return the number of equations
-    virtual size_t neq() {
+    virtual size_t neq() const {
         return m_nv;
     }
 
diff --git a/include/cantera/numerics/CVodesIntegrator.h b/include/cantera/numerics/CVodesIntegrator.h
index 28963483f4..92109da467 100644
--- a/include/cantera/numerics/CVodesIntegrator.h
+++ b/include/cantera/numerics/CVodesIntegrator.h
@@ -90,6 +90,10 @@ class CVodesIntegrator : public Integrator
 private:
     void sensInit(double t0, FuncEval& func);
 
+    //! Check whether a CVODES method indicated an error. If so, throw an exception
+    //! containing the method name and the error code stashed by the cvodes_err() function.
+    void checkError(long flag, const string& ctMethod, const string& cvodesMethod) const;
+
     size_t m_neq = 0;
     void* m_cvode_mem = nullptr;
     SundialsContext m_sundials_ctx; //!< SUNContext object for Sundials>=6.0
@@ -97,8 +101,16 @@ class CVodesIntegrator : public Integrator
     void* m_linsol_matrix = nullptr; //!< matrix used by Sundials
     FuncEval* m_func = nullptr;
     double m_t0 = 0.0;
-    double m_time; //!< The current integrator time
+
+    //! The current system time, corresponding to #m_y
+    double m_time;
+
+    //! The latest time reached by the integrator. May be greater than #m_time.
+    double m_tInteg;
+
+    //! The system state at #m_time
     N_Vector m_y = nullptr;
+
     N_Vector m_abstol = nullptr;
     N_Vector m_dky = nullptr;
     string m_type = "DENSE";
diff --git a/include/cantera/numerics/FuncEval.h b/include/cantera/numerics/FuncEval.h
index 01502fa5b1..02a7193001 100644
--- a/include/cantera/numerics/FuncEval.h
+++ b/include/cantera/numerics/FuncEval.h
@@ -150,10 +150,10 @@ class FuncEval
     }
 
     //! Number of equations.
-    virtual size_t neq()=0;
+    virtual size_t neq() const = 0;
 
     //! Number of sensitivity parameters.
-    virtual size_t nparams() {
+    virtual size_t nparams() const {
         return m_sens_params.size();
     }
 
diff --git a/include/cantera/numerics/IdasIntegrator.h b/include/cantera/numerics/IdasIntegrator.h
index 7096d9274b..5560dec0bf 100644
--- a/include/cantera/numerics/IdasIntegrator.h
+++ b/include/cantera/numerics/IdasIntegrator.h
@@ -107,7 +107,11 @@ class IdasIntegrator : public Integrator
     FuncEval* m_func = nullptr;
 
     double m_t0 = 0.0; //!< The start time for the integrator
-    double m_time; //!< The current integrator time
+    double m_time; //!< The current integrator time, corresponding to #m_y
+
+    //! The latest time reached by the integrator. May be greater than #m_time
+    double m_tInteg;
+
     N_Vector m_y = nullptr; //!< The current system state
     N_Vector m_ydot = nullptr; //!< The time derivatives of the system state
     N_Vector m_abstol = nullptr; //!< Absolute tolerances for each state variable
diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h
index 5068de5889..c4a3516b94 100644
--- a/include/cantera/zeroD/Reactor.h
+++ b/include/cantera/zeroD/Reactor.h
@@ -157,7 +157,7 @@ class Reactor : public ReactorBase
 
     //! Number of sensitivity parameters associated with this reactor
     //! (including walls)
-    virtual size_t nSensParams();
+    virtual size_t nSensParams() const;
 
     //! Add a sensitivity parameter associated with the reaction number *rxn*
     //! (in the homogeneous phase).
@@ -183,14 +183,14 @@ class Reactor : public ReactorBase
 
     //! Check whether Reactor object uses advance limits
     //! @returns           True if at least one limit is set, False otherwise
-    bool hasAdvanceLimits() {
+    bool hasAdvanceLimits() const {
         return !m_advancelimits.empty();
     }
 
     //! Retrieve absolute step size limits during advance
     //! @param[out] limits array of step size limits with length neq
     //! @returns           True if at least one limit is set, False otherwise
-    bool getAdvanceLimits(double* limits);
+    bool getAdvanceLimits(double* limits) const;
 
     //! Set individual step size limit for component name *nm*
     //! @param nm component name
diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h
index 9b21a54d53..369f7b82d3 100644
--- a/include/cantera/zeroD/ReactorNet.h
+++ b/include/cantera/zeroD/ReactorNet.h
@@ -60,7 +60,7 @@ class ReactorNet : public FuncEval
     }
 
     //! Get the maximum integrator step.
-    double maxTimeStep() {
+    double maxTimeStep() const {
         return m_maxstep;
     }
 
@@ -205,11 +205,11 @@ class ReactorNet : public FuncEval
                       double* ydot, double* p, Array2D* j);
 
     // overloaded methods of class FuncEval
-    virtual size_t neq() {
+    virtual size_t neq() const {
         return m_nv;
     }
 
-    size_t nReactors() {
+    size_t nReactors() const {
         return m_reactors.size();
     }
 
@@ -227,7 +227,7 @@ class ReactorNet : public FuncEval
 
     virtual void getConstraints(double* constraints);
 
-    virtual size_t nparams() {
+    virtual size_t nparams() const {
         return m_sens_params.size();
     }
 
@@ -254,7 +254,7 @@ class ReactorNet : public FuncEval
                                         double scale);
 
     //! The name of the p-th sensitivity parameter added to this ReactorNet.
-    const std::string& sensitivityParameterName(size_t p) {
+    const string& sensitivityParameterName(size_t p) const {
         return m_paramNames.at(p);
     }
 
@@ -284,10 +284,10 @@ class ReactorNet : public FuncEval
     void setAdvanceLimits(const double* limits);
 
     //! Check whether ReactorNet object uses advance limits
-    bool hasAdvanceLimits();
+    bool hasAdvanceLimits() const;
 
     //! Retrieve absolute step size limits during advance
-    bool getAdvanceLimits(double* limits);
+    bool getAdvanceLimits(double* limits) const;
 
     virtual void preconditionerSetup(double t, double* y, double gamma);
 
@@ -302,7 +302,7 @@ class ReactorNet : public FuncEval
 
 protected:
     //! Check that preconditioning is supported by all reactors in the network
-    virtual void checkPreconditionerSupported();
+    virtual void checkPreconditionerSupported() const;
 
     //! Update the preconditioner based on the already computed jacobian values
     virtual void updatePreconditioner(double gamma);
@@ -315,7 +315,7 @@ class ReactorNet : public FuncEval
     //! Returns the order used for last solution step of the ODE integrator
     //! The function is intended for internal use by ReactorNet::advance
     //! and deliberately not exposed in external interfaces.
-    virtual int lastOrder();
+    virtual int lastOrder() const;
 
     std::vector<Reactor*> m_reactors;
     std::unique_ptr<Integrator> m_integ;
diff --git a/samples/cxx/custom/custom.cpp b/samples/cxx/custom/custom.cpp
index e3ec0486ec..ca3784a476 100644
--- a/samples/cxx/custom/custom.cpp
+++ b/samples/cxx/custom/custom.cpp
@@ -113,7 +113,7 @@ class ReactorODEs : public FuncEval {
      * Number of equations in the ODE system.
      *   - overridden from FuncEval, called by the integrator during initialization.
      */
-    size_t neq() {
+    size_t neq() const {
         return m_nEqs;
     }
 
diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp
index eebaa76c3c..bb4aa6b301 100644
--- a/src/numerics/CVodesIntegrator.cpp
+++ b/src/numerics/CVodesIntegrator.cpp
@@ -236,10 +236,8 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func)
 
     int flag = CVodeSensInit(m_cvode_mem, static_cast<int>(m_np),
                              CV_STAGGERED, CVSensRhsFn(0), m_yS);
+    checkError(flag, "sensInit", "CVodeSensInit");
 
-    if (flag != CV_SUCCESS) {
-        throw CanteraError("CVodesIntegrator::sensInit", "Error in CVodeSensInit");
-    }
     vector_fp atol(m_np);
     for (size_t n = 0; n < m_np; n++) {
         // This scaling factor is tuned so that reaction and species enthalpy
@@ -247,6 +245,7 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func)
         atol[n] = m_abstolsens / func.m_paramScales[n];
     }
     flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data());
+    checkError(flag, "sensInit", "CVodeSensSStolerances");
 }
 
 void CVodesIntegrator::initialize(double t0, FuncEval& func)
@@ -254,6 +253,7 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
     m_neq = func.neq();
     m_t0 = t0;
     m_time = t0;
+    m_tInteg = t0;
     m_func = &func;
     func.clearErrors();
     // Initialize preconditioner if applied
@@ -310,39 +310,24 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
                                "CVodeInit failed.");
         }
     }
-    CVodeSetErrHandlerFn(m_cvode_mem, &cvodes_err, this);
+    flag = CVodeSetErrHandlerFn(m_cvode_mem, &cvodes_err, this);
 
     if (m_itol == CV_SV) {
         flag = CVodeSVtolerances(m_cvode_mem, m_reltol, m_abstol);
+        checkError(flag, "initialize", "CVodeSVtolerances");
     } else {
         flag = CVodeSStolerances(m_cvode_mem, m_reltol, m_abstols);
-    }
-    if (flag != CV_SUCCESS) {
-        if (flag == CV_MEM_FAIL) {
-            throw CanteraError("CVodesIntegrator::initialize",
-                               "Memory allocation failed.");
-        } else if (flag == CV_ILL_INPUT) {
-            throw CanteraError("CVodesIntegrator::initialize",
-                               "Illegal value for CVodeInit input argument.");
-        } else {
-            throw CanteraError("CVodesIntegrator::initialize",
-                               "CVodeInit failed.");
-        }
+        checkError(flag, "initialize", "CVodeSStolerances");
     }
 
     flag = CVodeSetUserData(m_cvode_mem, &func);
-    if (flag != CV_SUCCESS) {
-        throw CanteraError("CVodesIntegrator::initialize",
-                           "CVodeSetUserData failed.");
-    }
+    checkError(flag, "initialize", "CVodeSetUserData");
+
     if (func.nparams() > 0) {
         sensInit(t0, func);
         flag = CVodeSetSensParams(m_cvode_mem, func.m_sens_params.data(),
                                   func.m_paramScales.data(), NULL);
-        if (flag != CV_SUCCESS) {
-            throw CanteraError("CVodesIntegrator::initialize",
-                               "CVodeSetSensParams failed.");
-        }
+        checkError(flag, "initialize", "CVodeSetSensParams");
     }
     applyOptions();
 }
@@ -351,6 +336,7 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
 {
     m_t0 = t0;
     m_time = t0;
+    m_tInteg = t0;
     func.getState(NV_DATA_S(m_y));
     m_func = &func;
     func.clearErrors();
@@ -359,10 +345,7 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
         m_preconditioner->initialize(m_neq);
     }
     int result = CVodeReInit(m_cvode_mem, m_t0, m_y);
-    if (result != CV_SUCCESS) {
-        throw CanteraError("CVodesIntegrator::reinitialize",
-                           "CVodeReInit failed. result = {}", result);
-    }
+    checkError(result, "reinitialize", "CVodeReInit");
     applyOptions();
 }
 
@@ -492,7 +475,8 @@ void CVodesIntegrator::applyOptions()
     }
 
     if (m_maxord > 0) {
-        CVodeSetMaxOrd(m_cvode_mem, m_maxord);
+        int flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
+        checkError(flag, "applyOptions", "CVodeSetMaxOrd");
     }
     if (m_maxsteps > 0) {
         CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
@@ -512,25 +496,43 @@ void CVodesIntegrator::integrate(double tout)
 {
     if (tout == m_time) {
         return;
-    }
-    int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_NORMAL);
-    if (flag != CV_SUCCESS) {
-        string f_errs = m_func->getErrors();
-        if (!f_errs.empty()) {
-            f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
-        }
+    } else if (tout < m_time) {
         throw CanteraError("CVodesIntegrator::integrate",
-            "CVodes error encountered. Error code: {}\n{}\n"
-            "{}"
-            "Components with largest weighted error estimates:\n{}",
-            flag, m_error_message, f_errs, getErrorInfo(10));
+                           "Cannot integrate backwards in time.\n"
+                           "Requested time {} < current time {}",
+                           tout, m_time);
+    }
+    int nsteps = 0;
+    while (m_tInteg < tout) {
+        if (nsteps >= m_maxsteps) {
+            throw CanteraError("CVodesIntegrator::integrate",
+                "Maximum number of timesteps ({}) taken without reaching output "
+                "time ({}).\nCurrent integrator time: {}",
+                nsteps, tout, m_tInteg);
+        }
+        int flag = CVode(m_cvode_mem, tout, m_y, &m_tInteg, CV_ONE_STEP);
+        if (flag != CV_SUCCESS) {
+            string f_errs = m_func->getErrors();
+            if (!f_errs.empty()) {
+                f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
+            }
+            throw CanteraError("CVodesIntegrator::integrate",
+                "CVodes error encountered. Error code: {}\n{}\n"
+                "{}"
+                "Components with largest weighted error estimates:\n{}",
+                flag, m_error_message, f_errs, getErrorInfo(10));
+        }
+        nsteps++;
     }
+    int flag = CVodeGetDky(m_cvode_mem, tout, 0, m_y);
+    checkError(flag, "integrate", "CVodeGetDky");
+    m_time = tout;
     m_sens_ok = false;
 }
 
 double CVodesIntegrator::step(double tout)
 {
-    int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_ONE_STEP);
+    int flag = CVode(m_cvode_mem, tout, m_y, &m_tInteg, CV_ONE_STEP);
     if (flag != CV_SUCCESS) {
         string f_errs = m_func->getErrors();
         if (!f_errs.empty()) {
@@ -544,22 +546,14 @@ double CVodesIntegrator::step(double tout)
 
     }
     m_sens_ok = false;
+    m_time = m_tInteg;
     return m_time;
 }
 
 double* CVodesIntegrator::derivative(double tout, int n)
 {
     int flag = CVodeGetDky(m_cvode_mem, tout, n, m_dky);
-    if (flag != CV_SUCCESS) {
-        string f_errs = m_func->getErrors();
-        if (!f_errs.empty()) {
-            f_errs = "Exceptions caught evaluating derivative:\n" + f_errs;
-        }
-        throw CanteraError("CVodesIntegrator::derivative",
-            "CVodes error encountered. Error code: {}\n{}\n"
-            "{}",
-            flag, m_error_message, f_errs);
-    }
+    checkError(flag, "derivative", "CVodeGetDky");
     return NV_DATA_S(m_dky);
 }
 
@@ -643,11 +637,8 @@ double CVodesIntegrator::sensitivity(size_t k, size_t p)
         return 0.0;
     }
     if (!m_sens_ok && m_np) {
-        int flag = CVodeGetSens(m_cvode_mem, &m_time, m_yS);
-        if (flag != CV_SUCCESS) {
-            throw CanteraError("CVodesIntegrator::sensitivity",
-                               "CVodeGetSens failed. Error code: {}", flag);
-        }
+        int flag = CVodeGetSensDky(m_cvode_mem, m_time, 0, m_yS);
+        checkError(flag, "sensitivity", "CVodeGetSens");
         m_sens_ok = true;
     }
 
@@ -687,4 +678,20 @@ string CVodesIntegrator::getErrorInfo(int N)
     return to_string(s);
 }
 
+void CVodesIntegrator::checkError(long flag, const string& ctMethod,
+                                  const string& cvodesMethod) const
+{
+    if (flag == CV_SUCCESS) {
+        return;
+    } else if (flag == CV_MEM_NULL) {
+        throw CanteraError("CVodesIntegrator::" + ctMethod,
+                           "CVODES integrator is not initialized");
+    } else {
+        const char* flagname = CVodeGetReturnFlagName(flag);
+        throw CanteraError("CVodesIntegrator::" + ctMethod,
+            "{} returned error code {} ({}):\n{}",
+            cvodesMethod, flag, flagname, m_error_message);
+    }
+}
+
 }
diff --git a/src/numerics/IdasIntegrator.cpp b/src/numerics/IdasIntegrator.cpp
index cad7e7d055..f302d221f1 100644
--- a/src/numerics/IdasIntegrator.cpp
+++ b/src/numerics/IdasIntegrator.cpp
@@ -227,6 +227,7 @@ void IdasIntegrator::initialize(double t0, FuncEval& func)
     m_neq = func.neq();
     m_t0 = t0;
     m_time = t0;
+    m_tInteg = t0;
     m_func = &func;
     func.clearErrors();
 
@@ -317,6 +318,7 @@ void IdasIntegrator::reinitialize(double t0, FuncEval& func)
 {
     m_t0 = t0;
     m_time = t0;
+    m_tInteg = t0;
     func.getStateDae(NV_DATA_S(m_y), NV_DATA_S(m_ydot));
     m_func = &func;
     func.clearErrors();
@@ -451,24 +453,42 @@ void IdasIntegrator::integrate(double tout)
 {
     if (tout == m_time) {
         return;
-    }
-    int flag = IDASolve(m_ida_mem, tout, &m_time, m_y, m_ydot, IDA_NORMAL);
-    if (flag != IDA_SUCCESS) {
-        string f_errs = m_func->getErrors();
-        if (!f_errs.empty()) {
-            f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
-        }
+    } else if (tout < m_time) {
         throw CanteraError("IdasIntegrator::integrate",
-            "IDA error encountered. Error code: {}\n{}\n"
-            "{}"
-            "Components with largest weighted error estimates:\n{}",
-            flag, m_error_message, f_errs, getErrorInfo(10));
+                           "Cannot integrate backwards in time.\n"
+                           "Requested time {} < current time {}",
+                           tout, m_time);
+    }
+    int nsteps = 0;
+    while (m_tInteg < tout) {
+        if (nsteps >= m_maxsteps) {
+            throw CanteraError("IdasIntegrator::integrate",
+                "Maximum number of timesteps ({}) taken without reaching output "
+                "time ({}).\nCurrent integrator time: {}",
+                nsteps, tout, m_time);
+        }
+        int flag = IDASolve(m_ida_mem, tout, &m_tInteg, m_y, m_ydot, IDA_ONE_STEP);
+        if (flag != IDA_SUCCESS) {
+            string f_errs = m_func->getErrors();
+            if (!f_errs.empty()) {
+                f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
+            }
+            throw CanteraError("IdasIntegrator::integrate",
+                "IDA error encountered. Error code: {}\n{}\n"
+                "{}"
+                "Components with largest weighted error estimates:\n{}",
+                flag, m_error_message, f_errs, getErrorInfo(10));
+        }
+        nsteps++;
     }
+    int flag = IDAGetDky(m_ida_mem, tout, 0, m_y);
+    checkError(flag, "integrate", "IDAGetDky");
+    m_time = tout;
 }
 
 double IdasIntegrator::step(double tout)
 {
-    int flag = IDASolve(m_ida_mem, tout, &m_time, m_y, m_ydot, IDA_ONE_STEP);
+    int flag = IDASolve(m_ida_mem, tout, &m_tInteg, m_y, m_ydot, IDA_ONE_STEP);
     if (flag != IDA_SUCCESS) {
         string f_errs = m_func->getErrors();
         if (!f_errs.empty()) {
@@ -481,6 +501,7 @@ double IdasIntegrator::step(double tout)
             flag, f_errs, m_error_message, getErrorInfo(10));
 
     }
+    m_time = m_tInteg;
     return m_time;
 }
 
@@ -491,7 +512,7 @@ double IdasIntegrator::sensitivity(size_t k, size_t p)
         return 0.0;
     }
     if (!m_sens_ok && m_np) {
-        int flag = IDAGetSens(m_ida_mem, &m_time, m_yS);
+        int flag = IDAGetSensDky(m_ida_mem, m_time, 0, m_yS);
         checkError(flag, "sensitivity", "IDAGetSens");
         m_sens_ok = true;
     }
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index 9be81e215b..9e1c18cb08 100644
--- a/src/zeroD/Reactor.cpp
+++ b/src/zeroD/Reactor.cpp
@@ -109,7 +109,7 @@ void Reactor::initialize(double t0)
     m_work.resize(maxnt);
 }
 
-size_t Reactor::nSensParams()
+size_t Reactor::nSensParams() const
 {
     size_t ns = m_sensParams.size();
     for (auto& S : m_surfaces) {
@@ -549,7 +549,7 @@ void Reactor::setAdvanceLimits(const double *limits)
     }
 }
 
-bool Reactor::getAdvanceLimits(double *limits)
+bool Reactor::getAdvanceLimits(double *limits) const
 {
     bool has_limit = hasAdvanceLimits();
     if (has_limit) {
diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp
index c579e14cc4..e793bea5e8 100644
--- a/src/zeroD/ReactorNet.cpp
+++ b/src/zeroD/ReactorNet.cpp
@@ -272,7 +272,7 @@ void ReactorNet::getEstimate(double time, int k, double* yest)
     }
 }
 
-int ReactorNet::lastOrder()
+int ReactorNet::lastOrder() const
 {
     if (m_integ) {
         return m_integ->lastOrder();
@@ -425,7 +425,7 @@ void ReactorNet::setAdvanceLimits(const double *limits)
     }
 }
 
-bool ReactorNet::hasAdvanceLimits()
+bool ReactorNet::hasAdvanceLimits() const
 {
     bool has_limit = false;
     for (size_t n = 0; n < m_reactors.size(); n++) {
@@ -434,7 +434,7 @@ bool ReactorNet::hasAdvanceLimits()
     return has_limit;
 }
 
-bool ReactorNet::getAdvanceLimits(double *limits)
+bool ReactorNet::getAdvanceLimits(double *limits) const
 {
     bool has_limit = false;
     for (size_t n = 0; n < m_reactors.size(); n++) {
@@ -569,7 +569,7 @@ void ReactorNet::updatePreconditioner(double gamma)
     precon->updatePreconditioner();
 }
 
-void ReactorNet::checkPreconditionerSupported() {
+void ReactorNet::checkPreconditionerSupported() const {
     // check for non-mole-based reactors and throw an error otherwise
     for (auto reactor : m_reactors) {
         if (!reactor->preconditionerSupported()) {
diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py
index 4bb8d37f79..9f6da7aa6b 100644
--- a/test/python/test_reactor.py
+++ b/test/python/test_reactor.py
@@ -221,7 +221,7 @@ def test_maxsteps(self):
         self.net.max_time_step = max_step_size
         self.net.max_steps = max_steps
         with self.assertRaisesRegex(
-                ct.CanteraError, 'mxstep steps taken before reaching tout'):
+                ct.CanteraError, 'Maximum number of timesteps'):
             self.net.advance(1e-04)
         self.assertLessEqual(self.net.time, max_steps * max_step_size)
         self.assertEqual(self.net.max_steps, max_steps)
@@ -348,6 +348,12 @@ def test_advance_limits_invalid(self):
         with pytest.raises(ct.CanteraError, match="No component named 'spam'"):
             self.r1.set_advance_limit("spam", 0.1)
 
+    def test_advance_reverse(self):
+        self.make_reactors(n_reactors=1)
+        self.net.advance(0.1)
+        with pytest.raises(ct.CanteraError, match="backwards in time"):
+            self.net.advance(0.09)
+
     def test_heat_transfer2(self):
         # Result should be the same if (m * cp) / (U * A) is held constant
         self.make_reactors(T1=300, T2=1000)
@@ -1531,6 +1537,15 @@ def make_reactors(self, gas, surf):
         sim = ct.ReactorNet([r])
         return r, rsurf, sim
 
+    def test_advance_reverse(self):
+        surf, gas = self.import_phases()
+        gas.TPX = 1500, 4000, 'NH3:1.0, SiF4:0.4'
+        r, rsurf, sim = self.make_reactors(gas, surf)
+
+        sim.advance(0.1)
+        with pytest.raises(ct.CanteraError, match="backwards in time"):
+            sim.advance(0.09)
+
     def test_no_mass_flow_rate(self):
         surf, gas = self.import_phases()
         r = ct.FlowReactor(gas)
@@ -1605,7 +1620,7 @@ def test_max_steps(self):
         sim.max_steps = 13
         assert sim.max_steps == 13
 
-        with pytest.raises(ct.CanteraError, match="mxstep steps taken"):
+        with pytest.raises(ct.CanteraError, match="Maximum number of timesteps"):
             sim.advance(0.1)
 
         assert sim.solver_stats['steps'] == 13
@@ -1624,11 +1639,12 @@ def test_max_time_step(self):
         surf.TP = gas.TP
         r, rsurf, sim = self.make_reactors(gas, surf)
 
-        sim.max_time_step = 0.2
+        sim.max_time_step = 0.002
         sim.advance(0.01)
+        sim.step()
         x1 = sim.distance
         x2 = sim.step()
-        dx_limit = 0.05 * (x2-x1)
+        dx_limit = 0.1 * (x2-x1)
 
         sim.max_time_step = dx_limit
         assert sim.max_time_step == dx_limit
@@ -1746,7 +1762,7 @@ def test_initial_condition_tolerances(self):
         # With few steps allowed, won't be able to reach steady-state
         r.inlet_surface_max_error_failures = 10
         r.inlet_surface_max_steps = 200
-        with pytest.raises(ct.CanteraError, match="mxstep steps taken"):
+        with pytest.raises(ct.CanteraError, match="Maximum number of timesteps"):
             sim.initialize()
 
         # Relaxing the tolerances will allow the integrator to reach steady-state