diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd
index bbf11098b72..9f88c4d40e2 100644
--- a/interfaces/cython/cantera/reactor.pxd
+++ b/interfaces/cython/cantera/reactor.pxd
@@ -46,6 +46,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         size_t componentIndex(string&) except +translate_exception
         string componentName(size_t) except +translate_exception
         void getState(span[double]) except +translate_exception
+        void getStateDae(span[double], span[double]) except +translate_exception
         void setInitialVolume(double) except +translate_exception
         void addSensitivityReaction(size_t) except +translate_exception
 
@@ -190,7 +191,8 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         cbool verbose()
         void setVerbose(cbool)
         size_t neq()
-        void getState(span[double])
+        void getState(span[double]) except +translate_exception
+        void getStateDae(span[double], span[double]) except +translate_exception
         void getDerivative(int, span[double]) except +translate_exception
         void setAdvanceLimits(span[double])
         cbool getAdvanceLimits(span[double])
diff --git a/interfaces/cython/cantera/reactor.pyi b/interfaces/cython/cantera/reactor.pyi
index e175da66f35..b7dfd8259c5 100644
--- a/interfaces/cython/cantera/reactor.pyi
+++ b/interfaces/cython/cantera/reactor.pyi
@@ -43,6 +43,7 @@ class ReactorBase:
     @property
     def n_vars(self) -> int: ...
     def get_state(self) -> Array: ...
+    def get_state_dae(self) -> tuple[Array, Array]: ...
     def syncState(self) -> None: ...
     @property
     def phase(self) -> Solution: ...
@@ -157,7 +158,7 @@ class ExtensibleIdealGasMoleReactor(ExtensibleReactor): ...
 class ExtensibleConstPressureMoleReactor(ExtensibleReactor): ...
 class ExtensibleIdealGasConstPressureMoleReactor(ExtensibleReactor): ...
 
-class ReactorSurface:
+class ReactorSurface(ReactorBase):
     reactor_type: ClassVar[str]
     def __init__(
         self,
@@ -174,8 +175,6 @@ class ReactorSurface:
     @area.setter
     def area(self, A: float) -> None: ...
     @property
-    def phase(self) -> Solution: ...
-    @property
     def coverages(self) -> Array: ...
     @coverages.setter
     def coverages(self, coverages: Array) -> None: ...
@@ -472,6 +471,7 @@ class ReactorNet:
     @property
     def n_vars(self) -> int: ...
     def get_state(self) -> Array: ...
+    def get_state_dae(self) -> tuple[Array, Array]: ...
     def get_derivative(self, k: int) -> Array: ...
     @property
     def advance_limits(self) -> Array: ...
diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
index 20bd76bbcf5..6030b56307f 100644
--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -111,6 +111,19 @@ cdef class ReactorBase:
         self.rbase.getState(span[double](&y[0], y.size))
         return y
 
+    def get_state_dae(self):
+        """
+        Get the state vector and its time derivative for reactors formulated as DAEs
+        (namely, `FlowReactor`).
+
+        .. versionadded:: 4.0
+        """
+        cdef np.ndarray[np.double_t, ndim=1] y = np.zeros(self.n_vars)
+        cdef np.ndarray[np.double_t, ndim=1] yp = np.zeros(self.n_vars)
+        self.rbase.getStateDae(span[double](&y[0], y.size),
+                               span[double](&yp[0], yp.size))
+        return y, yp
+
     def syncState(self):
         """
         Set the state of the Reactor to match that of the associated
@@ -1966,6 +1979,22 @@ cdef class ReactorNet:
         self.net.getState(span[double](&y[0], y.size))
         return y
 
+    def get_state_dae(self):
+        """
+        Get the combined state vector and its time derivative for reactor networks
+        formulated as DAEs (namely, those containing `FlowReactor`).
+
+        The combined state vector consists of the concatenated state vectors of
+        all entities contained.
+
+        .. versionadded:: 4.0
+        """
+        cdef np.ndarray[np.double_t, ndim=1] y = np.zeros(self.n_vars)
+        cdef np.ndarray[np.double_t, ndim=1] yp = np.zeros(self.n_vars)
+        self.net.getStateDae(span[double](&y[0], y.size),
+                             span[double](&yp[0], yp.size))
+        return y, yp
+
     def get_derivative(self, k):
         """
         Get the k-th derivative of the state vector of the reactor network with respect
diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py
index 3b7ac7ba89d..04b6671d6f6 100644
--- a/test/python/test_reactor.py
+++ b/test/python/test_reactor.py
@@ -1904,6 +1904,10 @@ def test_reacting(self):
         r = ct.FlowReactor(g)
         r.mass_flow_rate = 10
         net = ct.ReactorNet([r])
+        y0, yp0 = net.get_state_dae()
+        assert y0[1] == approx(10 / g.density)
+        assert yp0[r.component_index('CO')] < 0
+        assert yp0[r.component_index('OH')] > 0
 
         i = 0
         while net.distance < 1.0:
@@ -1954,6 +1958,14 @@ def test_catalytic_surface(self):
         assert r.phase.density * r.area * r.speed == approx(mdot)
         assert sum(r.phase.Y) == approx(1.0)
 
+        y_surf, yp_surf = rsurf.get_state_dae()
+        y_gas, yp_gas = r.get_state_dae()
+        y_net, yp_net = sim.get_state_dae()
+        assert y_net[:r.n_vars] == approx(y_gas)
+        assert y_net[r.n_vars:] == approx(y_surf)
+        assert yp_net[:r.n_vars] == approx(yp_gas, abs=1e-4)
+        assert yp_net[r.n_vars:] == approx(yp_surf, abs=1e-4)
+
         sim.advance(1e-5)
         X = r.phase['CH4', 'H2', 'CO'].X
         assert X == approx([0.07748481, 0.048165072, 0.01446654])