diff --git a/include/cantera/numerics/CVodesIntegrator.h b/include/cantera/numerics/CVodesIntegrator.h
index 1c2c1d5c8..1332336f0 100644
--- a/include/cantera/numerics/CVodesIntegrator.h
+++ b/include/cantera/numerics/CVodesIntegrator.h
@@ -3,7 +3,7 @@
  */
 
 // This file is part of Cantera. See License.txt in the top-level directory or
-// at http://www.cantera.org/license.txt for license and copyright information.
+// at https://cantera.org/license.txt for license and copyright information.
 
 #ifndef CT_CVODESWRAPPER_H
 #define CT_CVODESWRAPPER_H
@@ -41,6 +41,8 @@ public:
     virtual doublereal step(double tout);
     virtual double& solution(size_t k);
     virtual double* solution();
+    virtual double* derivative(double tout, int n);
+    virtual int lastOrder() const;
     virtual int nEquations() const {
         return static_cast<int>(m_neq);
     }
@@ -89,6 +91,7 @@ private:
     double m_t0;
     double m_time; //!< The current integrator time
     N_Vector m_y, m_abstol;
+    N_Vector m_dky;
     int m_type;
     int m_itol;
     int m_method;
diff --git a/include/cantera/numerics/Integrator.h b/include/cantera/numerics/Integrator.h
index 6e0cdfc2c..272468ae9 100644
--- a/include/cantera/numerics/Integrator.h
+++ b/include/cantera/numerics/Integrator.h
@@ -7,7 +7,7 @@
  */
 
 // This file is part of Cantera. See License.txt in the top-level directory or
-// at http://www.cantera.org/license.txt for license and copyright information.
+// at https://cantera.org/license.txt for license and copyright information.
 
 #ifndef CT_INTEGRATOR_H
 #define CT_INTEGRATOR_H
@@ -141,6 +141,18 @@ public:
         return 0;
     }
 
+    //! n-th derivative of the output function at time tout.
+    virtual double* derivative(double tout, int n) {
+        warn("derivative");
+        return 0;
+    }
+
+    //! Order used during the last solution step
+    virtual int lastOrder() const {
+        warn("lastOrder");
+        return 0;
+    }
+
     //! The number of equations.
     virtual int nEquations() const {
         warn("nEquations");
diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h
index 168a13135..54ea826c0 100644
--- a/include/cantera/zeroD/Reactor.h
+++ b/include/cantera/zeroD/Reactor.h
@@ -132,6 +132,20 @@ public:
     //! @see componentIndex()
     virtual std::string componentName(size_t k);
 
+    //! Set absolute step size limits during advance
+    //! @param limits array of step size limits with length neq
+    virtual void setAdvanceLimits(const double* limits);
+
+    //! 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
+    virtual bool getAdvanceLimits(double* limits);
+
+    //! Set individual step size limit for compoment name *nm*
+    //! @param nm component name
+    //! @param limit value for step size limit
+    virtual void setAdvanceLimit(const std::string& nm, const double limit);
+
 protected:
     //! Set reaction rate multipliers based on the sensitivity variables in
     //! *params*.
@@ -180,6 +194,8 @@ protected:
     bool m_energy;
     size_t m_nv;
 
+    vector_fp m_advancelimits; //!< Advance step limit
+
     // Data associated each sensitivity parameter
     std::vector<SensitivityParameter> m_sensParams;
 };
diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h
index dedd8e1d2..4393f45c5 100644
--- a/include/cantera/zeroD/ReactorNet.h
+++ b/include/cantera/zeroD/ReactorNet.h
@@ -1,7 +1,7 @@
 //! @file ReactorNet.h
 
 // This file is part of Cantera. See License.txt in the top-level directory or
-// at http://www.cantera.org/license.txt for license and copyright information.
+// at https://cantera.org/license.txt for license and copyright information.
 
 #ifndef CT_REACTORNET_H
 #define CT_REACTORNET_H
@@ -81,6 +81,17 @@ public:
      */
     void advance(doublereal time);
 
+    /**
+     * Advance the state of all reactors in time. Take as many internal
+     * timesteps as necessary towards *time*. If *applylimit* is true,
+     * the advance step will be automatically reduced if needed to
+     * stay within limits (set by setAdvanceLimit).
+     * Returns the time at the end of integration.
+     * @param time Time to advance to (s).
+     * @param applylimit Limit advance step (boolean).
+     */
+    double advance(double time, bool applylimit);
+
     //! Advance the state of all reactors in time.
     double step();
 
@@ -164,6 +175,9 @@ public:
 
     virtual void getState(doublereal* y);
 
+    //! Return k-th derivative at the current time
+    virtual void getDerivative(int k, double* dky);
+
     virtual size_t nparams() {
         return m_sens_params.size();
     }
@@ -195,6 +209,10 @@ public:
         return m_paramNames.at(p);
     }
 
+    //! Initialize the reactor network. Called automatically the first time
+    //! advance or step is called.
+    void initialize();
+
     //! Reinitialize the integrator. Used to solve a new problem (different
     //! initial conditions) but with the same configuration of the reactor
     //! network. Can be called manually, or automatically after calling
@@ -222,10 +240,25 @@ public:
         return m_integ->maxSteps();
     }
 
+    //! Set absolute step size limits during advance
+    virtual void setAdvanceLimits(const double* limits);
+
+    //! Retrieve absolute step size limits during advance
+    virtual bool getAdvanceLimits(double* limits);
+
 protected:
-    //! Initialize the reactor network. Called automatically the first time
-    //! advance or step is called.
-    void initialize();
+
+    //! Estimate a future state based on current derivatives.
+    //! The function is intended for internal use by ReactorNet::advance
+    //! and deliberately not exposed in external interfaces.
+    virtual void getEstimate(double time, int k, double* yest);
+
+    //! 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() {
+        return m_integ->lastOrder();
+    }
 
     std::vector<Reactor*> m_reactors;
     std::unique_ptr<Integrator> m_integ;
@@ -251,6 +284,8 @@ protected:
     std::vector<std::string> m_paramNames;
 
     vector_fp m_ydot;
+    vector_fp m_yest;
+    vector_fp m_advancelimits;
 };
 }
 
diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd
index cdcc08bf0..d024f33df 100644
--- a/interfaces/cython/cantera/_cantera.pxd
+++ b/interfaces/cython/cantera/_cantera.pxd
@@ -524,7 +524,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         size_t neq()
         void getState(double*)
         void addSurface(CxxReactorSurface*)
-
+        void setAdvanceLimit(string&, double) except +translate_exception
         void addSensitivityReaction(size_t) except +translate_exception
         void addSensitivitySpeciesEnthalpy(size_t) except +translate_exception
         size_t nSensParams()
@@ -602,8 +602,9 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
     cdef cppclass CxxReactorNet "Cantera::ReactorNet":
         CxxReactorNet()
         void addReactor(CxxReactor&)
-        void advance(double) except +translate_exception
+        double advance(double, cbool) except +translate_exception
         double step() except +translate_exception
+        void initialize() except +translate_exception
         void reinitialize() except +translate_exception
         double time()
         void setInitialTime(double)
@@ -618,8 +619,11 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         void setVerbose(cbool)
         size_t neq()
         void getState(double*)
+        void getDerivative(int, double *) except +translate_exception
+        void setAdvanceLimits(double*)
+        cbool getAdvanceLimits(double*)
         string componentName(size_t) except +translate_exception
-
+        size_t globalComponentIndex(string&, int) except +translate_exception
         void setSensitivityTolerances(double, double)
         double rtolSensitivity()
         double atolSensitivity()
diff --git a/interfaces/cython/cantera/examples/reactors/reactor1.py b/interfaces/cython/cantera/examples/reactors/reactor1.py
index a7347a3dc..2ff2759ae 100644
--- a/interfaces/cython/cantera/examples/reactors/reactor1.py
+++ b/interfaces/cython/cantera/examples/reactors/reactor1.py
@@ -12,16 +12,22 @@ gas.TPX = 1001.0, ct.one_atm, 'H2:2,O2:1,N2:4'
 r = ct.IdealGasConstPressureReactor(gas)
 
 sim = ct.ReactorNet([r])
-time = 0.0
+sim.verbose = True
+
+# limit advance when temperature difference is exceeded
+delta_T_max = 20.
+r.set_advance_limit('temperature', delta_T_max)
+
+dt_max = 1.e-5
+t_end = 100 * dt_max
 states = ct.SolutionArray(gas, extra=['t'])
 
-print('%10s %10s %10s %14s' % ('t [s]','T [K]','P [Pa]','u [J/kg]'))
-for n in range(100):
-    time += 1.e-5
-    sim.advance(time)
-    states.append(r.thermo.state, t=time*1e3)
-    print('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T,
-                                           r.thermo.P, r.thermo.u))
+print('{:10s} {:10s} {:10s} {:14s}'.format('t [s]','T [K]','P [Pa]','u [J/kg]'))
+while sim.time < t_end:
+    sim.advance(sim.time + dt_max)
+    states.append(r.thermo.state, t=sim.time*1e3)
+    print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(sim.time, r.T,
+                                                       r.thermo.P, r.thermo.u))
 
 # Plot the results if matplotlib is installed.
 # See http://matplotlib.org/ to get it.
diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
index 92302359b..7d495271c 100644
--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -315,6 +315,17 @@ cdef class Reactor(ReactorBase):
         self.reactor.getState(&y[0])
         return y
 
+    def set_advance_limit(self, name, limit):
+        """
+        Limit absolute change of component *name* during `ReactorNet.advance`.
+        (positive *limit* values are considered; negative values disable a
+        previously set advance limit for a solution component). Note that
+        limits are disabled by default (with individual values set to -1.).
+        """
+        if limit is None:
+            limit = -1.
+        self.reactor.setAdvanceLimit(stringify(name), limit)
+
 
 cdef class Reservoir(ReactorBase):
     """
@@ -859,12 +870,17 @@ cdef class ReactorNet:
         self._reactors.append(r)
         self.net.addReactor(deref(r.reactor))
 
-    def advance(self, double t):
+    def advance(self, double t, pybool apply_limit=True):
         """
-        Advance the state of the reactor network in time from the current
-        time to time *t* [s], taking as many integrator timesteps as necessary.
+        Advance the state of the reactor network in time from the current time
+        towards time *t* [s], taking as many integrator timesteps as necessary.
+        If *apply_limit* is true and an advance limit is specified, the reactor
+        state at the end of the timestep is estimated prior to advancing. If
+        the difference exceed limits, the end time is reduced by half until
+        the projected end state remains within specified limits.
+        Returns the time reached at the end of integration.
         """
-        self.net.advance(t)
+        return self.net.advance(t, apply_limit)
 
     def step(self):
         """
@@ -873,6 +889,12 @@ cdef class ReactorNet:
         """
         return self.net.step()
 
+    def initialize(self):
+        """
+        Force initialization of the integrator after initial setup.
+        """
+        self.net.initialize()
+
     def reinitialize(self):
         """
         Reinitialize the integrator after making changing to the state of the
@@ -966,6 +988,17 @@ cdef class ReactorNet:
         def __set__(self, pybool v):
             self.net.setVerbose(v)
 
+    def global_component_index(self, name, int reactor):
+        """
+        Returns the index of a component named *name* of a reactor with index
+        *reactor* within the global state vector. I.e. this determines the
+        (absolute) index of the component, where *reactor* is the index of the
+        reactor that holds the component. *name* is either a species name or the
+        name of a reactor state variable, e.g. 'int_energy', 'temperature', etc.
+        depending on the reactor's equations.
+        """
+        return self.net.globalComponentIndex(stringify(name), reactor)
+
     def component_name(self, int i):
         """
         Return the name of the i-th component of the global state vector. The
@@ -1071,6 +1104,38 @@ cdef class ReactorNet:
         self.net.getState(&y[0])
         return y
 
+    def get_derivative(self, k):
+        """
+        Get the k-th time derivative of the state vector of the reactor network.
+        """
+        if not self.n_vars:
+            raise CanteraError('ReactorNet empty or not initialized.')
+        cdef np.ndarray[np.double_t, ndim = 1] dky = np.zeros(self.n_vars)
+        self.net.getDerivative(k, & dky[0])
+        return dky
+
+    property advance_limits:
+        """
+        Get or set absolute limits for state changes during `ReactorNet.advance`
+        (positive values are considered; negative values disable a previously
+        set advance limit for a solution component). Note that limits are
+        disabled by default (with individual values set to -1.).
+        """
+        def __get__(self):
+            cdef np.ndarray[np.double_t, ndim=1] limits = np.empty(self.n_vars)
+            self.net.getAdvanceLimits(&limits[0])
+            return limits
+
+        def __set__(self, limits):
+            if limits is None:
+                limits = -1. * np.ones([self.n_vars])
+            elif len(limits) != self.n_vars:
+                raise ValueError('array must be of length n_vars')
+
+            cdef np.ndarray[np.double_t, ndim=1] data = \
+                np.ascontiguousarray(limits, dtype=np.double)
+            self.net.setAdvanceLimits(&data[0])
+
     def advance_to_steady_state(self, int max_steps=10000,
                                 double residual_threshold=0., double atol=0.,
                                 pybool return_residuals=False):
diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py
index a7f7f201e..b39c77c3c 100644
--- a/interfaces/cython/cantera/test/test_reactor.py
+++ b/interfaces/cython/cantera/test/test_reactor.py
@@ -121,6 +121,22 @@ class TestReactor(utilities.CanteraTest):
         self.assertNear(P1, self.r1.thermo.P)
         self.assertNear(P2, self.r2.thermo.P)
 
+    def test_derivative(self):
+        T1, P1 = 300, 101325
+
+        self.make_reactors(n_reactors=1, T1=T1, P1=P1)
+        self.net.advance(1.0)
+
+        # compare cvode derivative to numerical derivative
+        dydt = self.net.get_derivative(1)
+        dt = -self.net.time
+        dy = -self.net.get_state()
+        self.net.step()
+        dt += self.net.time
+        dy += self.net.get_state()
+        for i in range(self.net.n_vars):
+            self.assertNear(dydt[i], dy[i]/dt)
+
     def test_timestepping(self):
         self.make_reactors()
 
@@ -205,6 +221,63 @@ class TestReactor(utilities.CanteraTest):
         self.assertTrue(n_baseline > n_rtol)
         self.assertTrue(n_baseline > n_atol)
 
+    def test_advance_limits(self):
+        P0 = 10 * ct.one_atm
+        T0 = 1100
+        X0 = 'H2:1.0, O2:0.5, AR:8.0'
+        self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0)
+
+        limit_H2 = .01
+        ix = self.net.global_component_index('H2', 0)
+        self.r1.set_advance_limit('H2', limit_H2)
+        self.assertEqual(self.net.advance_limits[ix], limit_H2)
+
+        self.r1.set_advance_limit('H2', None)
+        self.assertEqual(self.net.advance_limits[ix], -1.)
+
+        self.r1.set_advance_limit('H2', limit_H2)
+        self.net.advance_limits = None
+        self.assertEqual(self.net.advance_limits[ix], -1.)
+
+        self.r1.set_advance_limit('H2', limit_H2)
+        self.net.advance_limits = 0 * self.net.advance_limits - 1.
+        self.assertEqual(self.net.advance_limits[ix], -1.)
+
+    def test_advance_with_limits(self):
+        def integrate(limit_H2 = None, apply=True):
+            P0 = 10 * ct.one_atm
+            T0 = 1100
+            X0 = 'H2:1.0, O2:0.5, AR:8.0'
+            self.make_reactors(n_reactors=1, T1=T0, P1=P0, X1=X0)
+            if limit_H2 is not None:
+                self.r1.set_advance_limit('H2', limit_H2)
+                ix = self.net.global_component_index('H2', 0)
+                self.assertEqual(self.net.advance_limits[ix], limit_H2)
+
+            tEnd = 1.0
+            tStep = 1.e-3
+            nSteps = 0
+
+            t = tStep
+            while t < tEnd:
+                t_curr = self.net.advance(t, apply_limit=apply)
+                nSteps += 1
+                if t_curr == t:
+                    t += tStep
+
+            return nSteps
+
+        n_baseline = integrate()
+        n_advance_coarse = integrate(.01)
+        n_advance_fine = integrate(.001)
+        n_advance_negative = integrate(-.001)
+        n_advance_override = integrate(.001, False)
+
+        self.assertTrue(n_advance_coarse > n_baseline)
+        self.assertTrue(n_advance_fine > n_advance_coarse)
+        self.assertTrue(n_advance_negative == n_baseline)
+        self.assertTrue(n_advance_override == n_baseline)
+
     def test_heat_transfer1(self):
         # Connected reactors reach thermal equilibrium after some time
         self.make_reactors(T1=300, T2=1000)
diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp
index d1930aae6..806c67dd4 100644
--- a/src/numerics/CVodesIntegrator.cpp
+++ b/src/numerics/CVodesIntegrator.cpp
@@ -1,7 +1,7 @@
 //! @file CVodesIntegrator.cpp
 
 // This file is part of Cantera. See License.txt in the top-level directory or
-// at http://www.cantera.org/license.txt for license and copyright information.
+// at https://cantera.org/license.txt for license and copyright information.
 
 #include "cantera/numerics/CVodesIntegrator.h"
 #include "cantera/base/stringUtils.h"
@@ -85,6 +85,7 @@ CVodesIntegrator::CVodesIntegrator() :
     m_t0(0.0),
     m_y(0),
     m_abstol(0),
+    m_dky(0),
     m_type(DENSE+NOJAC),
     m_itol(CV_SS),
     m_method(CV_BDF),
@@ -126,6 +127,9 @@ CVodesIntegrator::~CVodesIntegrator()
     if (m_abstol) {
         N_VDestroy_Serial(m_abstol);
     }
+    if (m_dky) {
+        N_VDestroy_Serial(m_dky);
+    }
     if (m_yS) {
         N_VDestroyVectorArray_Serial(m_yS, static_cast<sd_size_t>(m_np));
     }
@@ -274,6 +278,11 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
     }
     m_y = N_VNew_Serial(static_cast<sd_size_t>(m_neq)); // allocate solution vector
     N_VConst(0.0, m_y);
+    if (m_dky) {
+        N_VDestroy_Serial(m_dky); // free derivative vector if already allocated
+    }
+    m_dky = N_VNew_Serial(static_cast<sd_size_t>(m_neq)); // allocate derivative vector
+    N_VConst(0.0, m_dky);
     // check abs tolerance array size
     if (m_itol == CV_SV && m_nabs < m_neq) {
         throw CanteraError("CVodesIntegrator::initialize",
@@ -484,6 +493,29 @@ double CVodesIntegrator::step(double tout)
     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);
+    }
+    return NV_DATA_S(m_dky);
+}
+
+int CVodesIntegrator::lastOrder() const
+{
+    int ord;
+    CVodeGetLastOrder(m_cvode_mem, &ord);
+    return ord;
+}
+
 int CVodesIntegrator::nEvals() const
 {
     long int ne;
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index 3ea29d7b9..924126de2 100644
--- a/src/zeroD/Reactor.cpp
+++ b/src/zeroD/Reactor.cpp
@@ -103,6 +103,7 @@ void Reactor::initialize(doublereal t0)
         }
     }
     m_work.resize(maxnt);
+    m_advancelimits.resize(m_nv, -1.0);
 }
 
 size_t Reactor::nSensParams()
@@ -429,4 +430,50 @@ void Reactor::resetSensitivity(double* params)
     }
 }
 
+void Reactor::setAdvanceLimits(const double *limits)
+{
+    if (m_thermo == 0) {
+        throw CanteraError("getState",
+                           "Error: reactor is empty.");
+    }
+    for (size_t j = 0; j < m_nv; j++) {
+        m_advancelimits[j] = limits[j];
+    }
+}
+
+bool Reactor::getAdvanceLimits(double *limits)
+{
+    bool has_limit = false;
+
+    for (size_t j = 0; j < m_nv; j++) {
+        limits[j] = m_advancelimits[j];
+        has_limit |= (limits[j] > 0.);
+    }
+
+    return has_limit;
+}
+
+void Reactor::setAdvanceLimit(const string& nm, const double limit)
+{
+    size_t k = componentIndex(nm);
+
+    if (m_thermo == 0) {
+        throw CanteraError("getState",
+                           "Error: reactor is empty.");
+    }
+    if (m_nv == 0) {
+        if (m_net == 0) {
+            throw CanteraError("Reactor::setAdvanceLimit",
+                               "Cannot set limit on a reactor that is not "
+                               "assigned to a ReactorNet object.");
+        } else {
+            m_net->initialize();
+        }
+    } else if (k > m_nv) {
+        throw CanteraError("Reactor::setAdvanceLimit",
+                           "Index out of bounds.");
+    }
+    m_advancelimits[k] = limit;
+}
+
 }
diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp
index 590d51acb..389da7112 100644
--- a/src/zeroD/ReactorNet.cpp
+++ b/src/zeroD/ReactorNet.cpp
@@ -1,7 +1,7 @@
 //! @file ReactorNet.cpp
 
 // This file is part of Cantera. See License.txt in the top-level directory or
-// at http://www.cantera.org/license.txt for license and copyright information.
+// at https://cantera.org/license.txt for license and copyright information.
 
 #include "cantera/zeroD/ReactorNet.h"
 #include "cantera/zeroD/FlowDevice.h"
@@ -98,6 +98,8 @@ void ReactorNet::initialize()
     }
 
     m_ydot.resize(m_nv,0.0);
+    m_yest.resize(m_nv,0.0);
+    m_advancelimits.resize(m_nv,-1.0);
     m_atol.resize(neq());
     fill(m_atol.begin(), m_atol.end(), m_atols);
     m_integ->setTolerances(m_rtol, neq(), m_atol.data());
@@ -136,6 +138,60 @@ void ReactorNet::advance(doublereal time)
     updateState(m_integ->solution());
 }
 
+double ReactorNet::advance(double time, bool applylimit)
+{
+    if (!m_init) {
+        initialize();
+    } else if (!m_integrator_init) {
+        reinitialize();
+    }
+
+    if (!applylimit) {
+        // take full step
+        advance(time);
+        return time;
+    }
+
+    bool limitadvance = getAdvanceLimits(m_advancelimits.data());
+    if (!limitadvance) {
+        // take full step
+        advance(time);
+        return time;
+    }
+
+    // ensure that gradient is available
+    while (lastOrder() < 1) {
+        step();
+    }
+
+    int k = lastOrder();
+    double t = time, delta;
+    double* y = m_integ->solution();
+
+    // reduce time step if limits are exceeded
+    while (true) {
+        bool exceeded = false;
+        getEstimate(t, k, &m_yest[0]);
+        for (size_t j = 0; j < m_nv; j++) {
+            delta = abs(m_yest[j] - y[j]);
+            if ( (m_advancelimits[j] > 0.) && ( delta > m_advancelimits[j]) ) {
+                exceeded = true;
+                if (m_verbose) {
+                    writelog("    Limiting global state vector component {:d} (dt = {:9.4g}):"
+                             "{:11.6g} > {:9.4g}\n",
+                             j, t - m_time, delta, m_advancelimits[j]);
+                }
+            }
+        }
+        if (!exceeded) {
+            break;
+        }
+        t = .5 * (m_time + t);
+    }
+    advance(t);
+    return t;
+}
+
 double ReactorNet::step()
 {
     if (!m_init) {
@@ -148,6 +204,26 @@ double ReactorNet::step()
     return m_time;
 }
 
+void ReactorNet::getEstimate(double time, int k, double* yest)
+{
+    // initialize
+    double* cvode_dky = m_integ->solution();
+    for (size_t j = 0; j < m_nv; j++) {
+        yest[j] = cvode_dky[j];
+    }
+
+    // Taylor expansion
+    double factor = 1.;
+    double deltat = time - m_time;
+    for (int n = 1; n <= k; n++) {
+        factor *= deltat / n;
+        cvode_dky = m_integ->derivative(m_time, n);
+        for (size_t j = 0; j < m_nv; j++) {
+            yest[j] += factor * cvode_dky[j];
+        }
+    }
+}
+
 void ReactorNet::addReactor(Reactor& r)
 {
     r.setNetwork(this);
@@ -218,6 +294,35 @@ void ReactorNet::getState(double* y)
     }
 }
 
+void ReactorNet::getDerivative(int k, double* dky)
+{
+    double* cvode_dky = m_integ->derivative(m_time, k);
+    for (size_t j = 0; j < m_nv; j++) {
+        dky[j] = cvode_dky[j];
+    }
+}
+
+void ReactorNet::setAdvanceLimits(const double *limits)
+{
+    if (!m_init) {
+        initialize();
+    }
+    for (size_t n = 0; n < m_reactors.size(); n++) {
+        m_reactors[n]->setAdvanceLimits(limits + m_start[n]);
+    }
+}
+
+bool ReactorNet::getAdvanceLimits(double *limits)
+{
+    bool has_limit = false;
+
+    for (size_t n = 0; n < m_reactors.size(); n++) {
+        has_limit |= m_reactors[n]->getAdvanceLimits(limits + m_start[n]);
+    }
+
+    return has_limit;
+}
+
 size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
 {
     if (!m_init) {