diff --git a/doc/sphinx/develop/index.md b/doc/sphinx/develop/index.md
index 8b19f52328..2cf4850171 100644
--- a/doc/sphinx/develop/index.md
+++ b/doc/sphinx/develop/index.md
@@ -1,6 +1,7 @@
 # Develop
 
 (sec-compiling)=
+
 ## Compiling Cantera from Source
 
 If you're interested in contributing new features to Cantera, or you want to try the
@@ -12,9 +13,9 @@ Cantera](compiling/configure-build) on your computer.
 
 The following additional references may also be useful:
 
-- [](compiling/dependencies.md)
-- [](compiling/config-options)
-- [](compiling/special-cases)
+-   [](compiling/dependencies.md)
+-   [](compiling/config-options)
+-   [](compiling/special-cases)
 
 ```{toctree}
 :caption: Compiling Cantera from Source
@@ -35,7 +36,7 @@ compiling/special-cases
 This section is a work in progress.
 ```
 
-- [](reactor-integration)
+-   [](reactor-integration)
 
 ```{toctree}
 :caption: How Cantera Works
diff --git a/include/cantera/base/Delegator.h b/include/cantera/base/Delegator.h
index 3a470e5e11..12122f5193 100644
--- a/include/cantera/base/Delegator.h
+++ b/include/cantera/base/Delegator.h
@@ -10,6 +10,7 @@
 #include "cantera/base/Units.h"
 #include "cantera/base/ctexceptions.h"
 #include "cantera/base/ExtensionManager.h"
+#include "cantera/numerics/eigen_sparse.h"
 #include <array>
 #include <list>
 
@@ -226,6 +227,20 @@ class Delegator
         *m_funcs_v_d_dp_dp[name] = makeDelegate(func, when, *m_funcs_v_d_dp_dp[name]);
     }
 
+    //! Set delegates for member functions with the signature
+    //! `void(vector<Eigen::Triplet<double>>&)`
+    void setDelegate(const string& name,
+                     const function<void(vector<Eigen::Triplet<double>>&)>& func,
+                     const string& when)
+    {
+        if (!m_funcs_v_vet.count(name)) {
+            throw NotImplementedError("Delegator::setDelegate",
+                "for function '{}' with signature "
+                "'void(vector<Eigen::Triplet<double>>&)'.", name);
+        }
+        *m_funcs_v_vet[name] = makeDelegate(func, when, *m_funcs_v_vet[name]);
+    }
+
     //! Set delegates for member functions with the signature
     //! `void(double*, double*, double*)`
     void setDelegate(
@@ -386,6 +401,16 @@ class Delegator
         m_funcs_v_dp_dp_dp[name] = &target;
     }
 
+    //! Install a function with the signature `void(vector<Eigen::Triplet<double>>&)`
+    //! as being delegatable
+    void install(const string& name,
+                 function<void(vector<Eigen::Triplet<double>>&)>& target,
+                 const function<void(vector<Eigen::Triplet<double>>&)>& base)
+    {
+        target = base;
+        m_funcs_v_vet[name] = &target;
+    }
+
     //! Install a function with the signature `double(void*)` as being delegatable
     void install(const string& name, function<double(void*)>& target,
                  const function<double(void*)>& func)
@@ -514,6 +539,7 @@ class Delegator
     //! - `sz` for `size_t`
     //! - `AM` for `AnyMap`
     //! - `US` for `UnitStack`
+    //! - `VET` for `vector<Eigen::Triplet<double>>`
     //! - prefix `c` for `const` arguments
     //! - suffix `r` for reference arguments
     //! - suffix `p` for pointer arguments
@@ -532,6 +558,7 @@ class Delegator
         function<void(std::array<size_t, 2>, double, double*, double*)>*> m_funcs_v_d_dp_dp;
     map<string,
         function<void(std::array<size_t, 3>, double*, double*, double*)>*> m_funcs_v_dp_dp_dp;
+    map<string, function<void(vector<Eigen::Triplet<double>>&)>*> m_funcs_v_vet;
 
     // Delegates with a return value
     map<string, function<double(void*)>> m_base_d_vp;
@@ -542,6 +569,7 @@ class Delegator
 
     map<string, function<size_t(const string&)>> m_base_sz_csr;
     map<string, function<size_t(const string&)>*> m_funcs_sz_csr;
+
     //! @}
 
     //! Handles to wrappers for the delegated object in external language interfaces.
diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h
index 34582a50b5..49cbe88f5a 100644
--- a/include/cantera/zeroD/FlowDevice.h
+++ b/include/cantera/zeroD/FlowDevice.h
@@ -10,6 +10,7 @@
 #include "cantera/base/global.h"
 #include "cantera/base/ctexceptions.h"
 #include "ConnectorNode.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -163,7 +164,42 @@ class FlowDevice : public ConnectorNode
         m_time = time;
     }
 
+    //! Build the Jacobian terms specific to the flow device for the given connected
+    //! reactor.
+    //! @param r a pointer to the calling reactor
+    //! @param jacVector a vector of triplets to be added to the jacobian for the
+    //! reactor
+    //! @warning This function is an experimental part of the %Cantera API and may be
+    //! changed
+    //! or removed without notice.
+    //! @since New in %Cantera 3.1.
+    //!
+    virtual void buildReactorJacobian(ReactorBase* r,
+        vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError(type() + "::buildReactorJacobian");
+    }
+
+    //! Build the Jacobian terms specific to the flow device for the network. These
+    //! terms
+    //! will be adjusted to the networks indexing system outside of the reactor.
+    //! @param jacVector a vector of triplets to be added to the jacobian for the
+    //! reactor
+    //! @warning This function is an experimental part of the %Cantera API and may be
+    //! changed
+    //! or removed without notice.
+    //! @since New in %Cantera 3.1.
+    //!
+    virtual void buildNetworkJacobian(vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError(type() + "::buildNetworkJacobian");
+    }
+
 protected:
+    string m_name;  //!< Flow device name.
+    bool m_defaultNameSet = false;  //!< `true` if default name has been previously set.
+    //! a variable to switch on and off so calculations are not doubled by the calling
+    //! reactor or network
+    bool m_jac_calculated = false;
+
     double m_mdot = Undef;
 
     //! Function set by setPressureFunction; used by updateMassFlowRate
diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h
index 27de6d57c4..62f75206a3 100644
--- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h
+++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h
@@ -39,10 +39,14 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor
     //! Neglects derivatives with respect to mole fractions that would generate a
     //! fully-dense Jacobian. Currently also neglects terms related to interactions
     //! between reactors, for example via inlets and outlets.
-    Eigen::SparseMatrix<double> jacobian() override;
+    void buildJacobian(vector<Eigen::Triplet<double>>& jacVector) override;
 
     bool preconditionerSupported() const override { return true; };
 
+    double temperature_ddni(size_t index) override;
+
+    size_t speciesOffset() const override { return m_sidx; };
+
     size_t componentIndex(const string& nm) const override;
     string componentName(size_t k) override;
     double upperBound(size_t k) const override;
diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h
index 9a2fceece4..48b6a0ccd6 100644
--- a/include/cantera/zeroD/IdealGasMoleReactor.h
+++ b/include/cantera/zeroD/IdealGasMoleReactor.h
@@ -40,15 +40,15 @@ class IdealGasMoleReactor : public MoleReactor
 
     void updateState(double* y) override;
 
-    //! Calculate an approximate Jacobian to accelerate preconditioned solvers
-
-    //! Neglects derivatives with respect to mole fractions that would generate a
-    //! fully-dense Jacobian. Currently, also neglects terms related to interactions
-    //! between reactors, for example via inlets and outlets.
-    Eigen::SparseMatrix<double> jacobian() override;
+    //! Calculate the Jacobian to accelerate preconditioned solvers
+    void buildJacobian(vector<Eigen::Triplet<double>>& jacVector) override;
 
     bool preconditionerSupported() const override {return true;};
 
+    double temperature_ddni(size_t index) override;
+
+    size_t speciesOffset() const override { return m_sidx; };
+
 protected:
     vector<double> m_uk; //!< Species molar internal energies
 };
diff --git a/include/cantera/zeroD/MoleReactor.h b/include/cantera/zeroD/MoleReactor.h
index e15e07695d..723d4b29df 100644
--- a/include/cantera/zeroD/MoleReactor.h
+++ b/include/cantera/zeroD/MoleReactor.h
@@ -7,6 +7,7 @@
 #define CT_MOLEREACTOR_H
 
 #include "Reactor.h"
+#include "cantera/zeroD/Wall.h"
 
 namespace Cantera
 {
@@ -40,6 +41,8 @@ class MoleReactor : public Reactor
     double lowerBound(size_t k) const override;
     void resetBadValues(double* y) override;
 
+    size_t energyIndex() const override { return m_eidx; };
+
 protected:
     //! For each surface in the reactor, update vector of triplets with all relevant
     //! surface jacobian derivatives of species with respect to species
@@ -62,6 +65,9 @@ class MoleReactor : public Reactor
 
     //! const value for the species start index
     const size_t m_sidx = 2;
+
+    //! index of state variable associated with energy
+    const size_t m_eidx = 0;
 };
 
 }
diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h
index fc5e48e976..84522aad79 100644
--- a/include/cantera/zeroD/Reactor.h
+++ b/include/cantera/zeroD/Reactor.h
@@ -77,12 +77,6 @@ class Reactor : public ReactorBase
         m_chem = cflag;
     }
 
-    //! Returns `true` if changes in the reactor composition due to chemical reactions
-    //! are enabled.
-    bool chemistryEnabled() const {
-        return m_chem;
-    }
-
     void setEnergy(int eflag=1) override {
         if (eflag > 0) {
             m_energy = true;
@@ -91,11 +85,6 @@ class Reactor : public ReactorBase
         }
     }
 
-    //! Returns `true` if solution of the energy equation is enabled.
-    bool energyEnabled() const {
-        return m_energy;
-    }
-
     //! Number of equations (state variables) for this reactor
     size_t neq() {
         if (!m_nv) {
@@ -209,7 +198,18 @@ class Reactor : public ReactorBase
     //! @param limit value for step size limit
     void setAdvanceLimit(const string& nm, const double limit);
 
+    //! A wrapper for the Jacobian function to return the Eigen::SparseMatrix<double>
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //! @ingroup derivGroup
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual Eigen::SparseMatrix<double> jacobian();
+
     //! Calculate the Jacobian of a specific Reactor specialization.
+    //! @param jacVector vector where jacobian triplets are added
     //! @warning Depending on the particular implementation, this may return an
     //! approximate Jacobian intended only for use in forming a preconditioner for
     //! iterative solvers.
@@ -217,10 +217,32 @@ class Reactor : public ReactorBase
     //!
     //! @warning  This method is an experimental part of the %Cantera
     //! API and may be changed or removed without notice.
-    virtual Eigen::SparseMatrix<double> jacobian() {
-        throw NotImplementedError("Reactor::jacobian");
+    virtual void buildJacobian(vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError(type() + "::buildJacobian");
     }
 
+    //! Calculate the Jacobian of a Reactor specialization for wall contributions.
+    //! @param jacVector vector where jacobian triplets are added
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //! @ingroup derivGroup
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual void buildWallJacobian(vector<Eigen::Triplet<double>>& jacVector);
+
+    //! Calculate flow contributions to the Jacobian of a Reactor specialization.
+    //! @param jacVector vector where jacobian triplets are added
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //! @ingroup derivGroup
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual void buildFlowJacobian(vector<Eigen::Triplet<double>>& jacVector);
+
     //! Calculate the reactor-specific Jacobian using a finite difference method.
     //!
     //! This method is used only for informational purposes. Jacobian calculations
@@ -303,8 +325,6 @@ class Reactor : public ReactorBase
 
     vector<double> m_wdot; //!< Species net molar production rates
     vector<double> m_uk; //!< Species molar internal energies
-    bool m_chem = false;
-    bool m_energy = true;
     size_t m_nv = 0;
     size_t m_nv_surf; //!!< Number of variables associated with reactor surfaces
 
@@ -312,6 +332,10 @@ class Reactor : public ReactorBase
 
     //! Vector of triplets representing the jacobian
     vector<Eigen::Triplet<double>> m_jac_trips;
+    //! Boolean to skip walls in jacobian
+    bool m_jac_skip_walls = false;
+    //! Boolean to skip flow devices in jacobian
+    bool m_jac_skip_flow_devices = false;
 };
 }
 
diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h
index 828d0bb68b..8601cb4529 100644
--- a/include/cantera/zeroD/ReactorBase.h
+++ b/include/cantera/zeroD/ReactorBase.h
@@ -320,6 +320,34 @@ class ReactorBase
         return m_sensParams.size();
     }
 
+    //! Calculate the derivative of T with respect to the ith species in the energy
+    //! conservation equation based on the reactor specific equation of state.
+    //! @param index index of the species the derivative is with respect too
+    //! @warning This function is an experimental part of the %Cantera API and may
+    //! be changed or removed without notice.
+    //! @since New in %Cantera 3.1.
+    //!
+    virtual double temperature_ddni(size_t index) {
+        throw NotImplementedError("Reactor::temperature_ddni");
+    }
+
+    //! Return the index associated with energy of the system
+    virtual size_t energyIndex() const { return m_eidx; };
+
+    //! Return the offset between species and state variables
+    virtual size_t speciesOffset() const { return m_sidx; };
+
+    //! Returns `true` if solution of the energy equation is enabled.
+    virtual bool energyEnabled() const {
+        return m_energy;
+    }
+
+    //! Returns `true` if changes in the reactor composition due to chemical reactions
+    //! are enabled.
+    bool chemistryEnabled() const {
+        return m_chem;
+    }
+
 protected:
     //! Specify the mixture contained in the reactor. Note that a pointer to
     //! this substance is stored, and as the integration proceeds, the state of
@@ -340,6 +368,12 @@ class ReactorBase
     //! Number of homogeneous species in the mixture
     size_t m_nsp = 0;
 
+    //! species offset in the state vector
+    const size_t m_sidx = 3;
+
+    //! index of state variable associated with energy
+    const size_t m_eidx = 1;
+
     ThermoPhase* m_thermo = nullptr;
     double m_vol = 0.0; //!< Current volume of the reactor [m^3]
     double m_mass = 0.0; //!< Current mass of the reactor [kg]
@@ -369,6 +403,12 @@ class ReactorBase
 
     // Data associated each sensitivity parameter
     vector<SensitivityParameter> m_sensParams;
+
+    //! A bool that enables the energy equation
+    bool m_energy = true;
+
+    //! A bool that enables the chemical kinetics equations
+    bool m_chem = false;
 };
 }
 
diff --git a/include/cantera/zeroD/ReactorDelegator.h b/include/cantera/zeroD/ReactorDelegator.h
index f2556a7871..5860e6a799 100644
--- a/include/cantera/zeroD/ReactorDelegator.h
+++ b/include/cantera/zeroD/ReactorDelegator.h
@@ -10,6 +10,7 @@
 #include "cantera/base/Delegator.h"
 #include "cantera/zeroD/ReactorSurface.h"
 #include "cantera/thermo/SurfPhase.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -51,6 +52,10 @@ class ReactorAccessor
     //! Set the state of the thermo object for surface *n* to correspond to the
     //! state of that surface
     virtual void restoreSurfaceState(size_t n) = 0;
+
+    //! Public access to the default evaluation function so it can be used in
+    //! replace functions
+    virtual void defaultEval(double t, double* LHS, double* RHS) = 0;
 };
 
 //! Delegate methods of the Reactor class to external functions
@@ -69,7 +74,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
         install("updateState", m_updateState,
             [this](std::array<size_t, 1> sizes, double* y) { R::updateState(y); });
         install("updateSurfaceState", m_updateSurfaceState,
-            [this](std::array<size_t, 1> sizes, double* y) { R::updateSurfaceState(y); });
+            [this](std::array<size_t, 1> sizes, double* y)
+            { R::updateSurfaceState(y); });
         install("getSurfaceInitialConditions", m_getSurfaceInitialConditions,
             [this](std::array<size_t, 1> sizes, double* y) {
                 R::getSurfaceInitialConditions(y);
@@ -84,8 +90,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
         );
         install("evalWalls", m_evalWalls, [this](double t) { R::evalWalls(t); });
         install("evalSurfaces", m_evalSurfaces,
-            [this](std::array<size_t, 3> sizes, double* LHS, double* RHS, double* sdot) {
-                R::evalSurfaces(LHS, RHS, sdot);
+            [this](std::array<size_t, 3> sizes, double* LHS, double* RHS, double* sd) {
+                R::evalSurfaces(LHS, RHS, sd);
             }
         );
         install("componentName", m_componentName,
@@ -94,6 +100,8 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
             [this](const string& nm) { return R::componentIndex(nm); });
         install("speciesIndex", m_speciesIndex,
             [this](const string& nm) { return R::speciesIndex(nm); });
+        install("buildJacobian", m_build_jacobian,
+            [this](vector<Eigen::Triplet<double>>& jv) { R::buildJacobian(jv); });
     }
 
     // Overrides of Reactor methods
@@ -160,6 +168,14 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
         return m_speciesIndex(nm);
     }
 
+    void buildJacobian(vector<Eigen::Triplet<double>>& jacVector) override {
+        m_build_jacobian(jacVector);
+    }
+
+    void defaultEval(double t, double* LHS, double* RHS) override {
+        R::eval(t, LHS, RHS);
+    }
+
     // Public access to protected Reactor variables needed by derived classes
 
     void setNEq(size_t n) override {
@@ -204,6 +220,7 @@ class ReactorDelegator : public Delegator, public R, public ReactorAccessor
     function<string(size_t)> m_componentName;
     function<size_t(const string&)> m_componentIndex;
     function<size_t(const string&)> m_speciesIndex;
+    function<void(vector<Eigen::Triplet<double>>&)> m_build_jacobian;
 };
 
 }
diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h
index cd2958341b..b00df70b98 100644
--- a/include/cantera/zeroD/ReactorNet.h
+++ b/include/cantera/zeroD/ReactorNet.h
@@ -7,6 +7,8 @@
 #define CT_REACTORNET_H
 
 #include "Reactor.h"
+#include "Wall.h"
+#include "FlowDevice.h"
 #include "cantera/numerics/FuncEval.h"
 #include "cantera/numerics/SteadyStateSystem.h"
 
@@ -286,6 +288,11 @@ class ReactorNet : public FuncEval
     //! reactor network.
     size_t globalComponentIndex(const string& component, size_t reactor=0);
 
+    //! Return the index corresponding to the start of the reactor specific state
+    //! vector in the reactor with index *reactor* in the global state vector for the
+    //! reactor network.
+    size_t globalStartIndex(ReactorBase* curr_reactor);
+
     //! Return the name of the i-th component of the global state vector. The
     //! name returned includes both the name of the reactor and the specific
     //! component, for example `'reactor1: CH4'`.
@@ -363,11 +370,48 @@ class ReactorNet : public FuncEval
     //! @param settings the settings map propagated to all reactors and kinetics objects
     virtual void setDerivativeSettings(AnyMap& settings);
 
+    //! Calculate the system Jacobian using a finite difference method.
+    //!
+    //! This method is used only for informational purposes. Jacobian calculations
+    //! for the full reactor system are handled internally by CVODES.
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual Eigen::SparseMatrix<double> finiteDifferenceJacobian();
+
+    //! A wrapper method to calculate the system jacobian as Eigen::SparseMatrix<double>
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual Eigen::SparseMatrix<double> jacobian() {
+        vector<Eigen::Triplet<double>> jac_trips;
+        // Add before, during, after evals
+        buildJacobian(jac_trips);
+        // construct jacobian from vector
+        Eigen::SparseMatrix<double> jac(m_nv, m_nv);
+        jac.setFromTriplets(jac_trips.begin(), jac_trips.end());
+        return jac;
+    }
+
 protected:
     //! Add the reactor *r* to this reactor network.
     //! @since  Changed in %Cantera 3.2. Previous version used a reference.
     void addReactor(shared_ptr<ReactorBase> reactor);
 
+    //! Calculate the Jacobian of the entire reactor network.
+    //! @param jacVector vector where jacobian triplets are added
+    //! @warning Depending on the particular implementation, this may return an
+    //! approximate Jacobian intended only for use in forming a preconditioner for
+    //! iterative solvers.
+    //! @ingroup derivGroup
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    virtual void buildJacobian(vector<Eigen::Triplet<double>>& jacVector);
+
     //! Check that preconditioning is supported by all reactors in the network
     virtual void checkPreconditionerSupported() const;
 
@@ -430,6 +474,14 @@ class ReactorNet : public FuncEval
     //! "left hand side" of each governing equation
     vector<double> m_LHS;
     vector<double> m_RHS;
+
+    //! derivative settings
+    bool m_jac_skip_walls = false;
+    bool m_jac_skip_flow_devices = false;
+    //! set to store walls for Jacobian calculation
+    set<WallBase*> m_walls;
+    //! set to store flow devices for Jacobian calculation
+    set<FlowDevice*> m_flow_devices;
 };
 
 
diff --git a/include/cantera/zeroD/Wall.h b/include/cantera/zeroD/Wall.h
index e88b3165a3..ee65eebf7e 100644
--- a/include/cantera/zeroD/Wall.h
+++ b/include/cantera/zeroD/Wall.h
@@ -9,6 +9,7 @@
 #include "cantera/base/ctexceptions.h"
 #include "cantera/zeroD/ReactorBase.h"
 #include "ConnectorNode.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -90,6 +91,30 @@ class WallBase : public ConnectorNode
         m_time = time;
     }
 
+    //! Build the Jacobian terms specific to the flow device for the given connected
+    //! reactor.
+    //! @param r a pointer to the calling reactor
+    //! @param jacVector a vector of triplets to be added to the reactor Jacobian
+    //! @warning This function is an experimental part of the %Cantera API and may be
+    //! changed or removed without notice.
+    //! @since New in %Cantera 3.1.
+    //!
+    virtual void buildReactorJacobian(ReactorBase* r,
+        vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError("WallBase::buildReactorJacobian");
+    }
+
+    //! Build the Jacobian cross-reactor terms specific to the flow device for the
+    //! network.
+    //! @param jacVector a vector of triplets to be added to the network Jacobian
+    //! @warning This function is an experimental part of the %Cantera API and may be
+    //! changed or removed without notice.
+    //! @since New in %Cantera 3.1.
+    //!
+    virtual void buildNetworkJacobian(vector<Eigen::Triplet<double>>& jacVector) {
+        throw NotImplementedError("WallBase::buildNetworkJacobian");
+    }
+
 protected:
     ReactorBase* m_left = nullptr;
     ReactorBase* m_right = nullptr;
@@ -227,6 +252,12 @@ class Wall : public WallBase
         return m_k;
     }
 
+    void buildReactorJacobian(ReactorBase* r,
+        vector<Eigen::Triplet<double>>& jacVector) override;
+
+    void buildNetworkJacobian(vector<Eigen::Triplet<double>>& jacVector)
+        override;
+
 protected:
 
     //! expansion rate coefficient
diff --git a/interfaces/cython/cantera/_utils.pxd b/interfaces/cython/cantera/_utils.pxd
index ade43f1e14..3949a296d4 100644
--- a/interfaces/cython/cantera/_utils.pxd
+++ b/interfaces/cython/cantera/_utils.pxd
@@ -9,6 +9,7 @@ from libcpp.memory cimport unique_ptr, make_unique
 
 from .ctcxx cimport *
 from .units cimport UnitSystem, CxxUnits
+from .delegator cimport CxxEigenTriplet
 
 cdef extern from "cantera/base/AnyMap.h" namespace "Cantera":
     cdef cppclass CxxAnyValue "Cantera::AnyValue"
@@ -118,3 +119,5 @@ cdef anymap_to_py(CxxAnyMap& m)
 
 cdef CxxAnyValue python_to_anyvalue(item, name=*) except *
 cdef anyvalue_to_python(string name, CxxAnyValue& v)
+
+cdef CxxEigenTriplet get_triplet(row, col, val) except *
diff --git a/interfaces/cython/cantera/_utils.pyx b/interfaces/cython/cantera/_utils.pyx
index aec6f25a0f..db2ece1694 100644
--- a/interfaces/cython/cantera/_utils.pyx
+++ b/interfaces/cython/cantera/_utils.pyx
@@ -5,6 +5,7 @@ import os
 import warnings
 from cpython.ref cimport PyObject
 from libcpp.utility cimport move
+from cython.operator cimport dereference
 import numbers as _numbers
 import importlib as _importlib
 from collections import namedtuple as _namedtuple
@@ -537,6 +538,10 @@ cdef vector[vector[string]] list2_string_to_anyvalue(data):
             v[i][j] = stringify(jtem)
     return v
 
+cdef CxxEigenTriplet get_triplet(row, col, val):
+    cdef CxxEigenTriplet* trip_ptr = new CxxEigenTriplet(row, col, val)
+    return dereference(trip_ptr)
+
 def _py_to_any_to_py(dd):
     # used for internal testing purposes only
     cdef string name = stringify("test")
diff --git a/interfaces/cython/cantera/delegator.pxd b/interfaces/cython/cantera/delegator.pxd
index 721ddd736a..09843a3faa 100644
--- a/interfaces/cython/cantera/delegator.pxd
+++ b/interfaces/cython/cantera/delegator.pxd
@@ -7,6 +7,15 @@
 from .ctcxx cimport *
 from .func1 cimport *
 from .units cimport CxxUnitStack
+# from .kinetics cimport *
+
+cdef extern from "cantera/numerics/eigen_sparse.h" namespace "Eigen":
+    cdef cppclass CxxEigenTriplet "Eigen::Triplet<double>":
+        CxxEigenTriplet()
+        CxxEigenTriplet(size_t, size_t, double)
+        size_t row()
+        size_t col()
+        size_t value()
 
 cdef extern from "<array>" namespace "std" nogil:
     cdef cppclass size_array1 "std::array<size_t, 1>":
@@ -52,7 +61,7 @@ cdef extern from "cantera/base/Delegator.h" namespace "Cantera":
         void setDelegate(string&, function[int(double&, void*)], string&) except +translate_exception
         void setDelegate(string&, function[int(string&, size_t)], string&) except +translate_exception
         void setDelegate(string&, function[int(size_t&, string&)], string&) except +translate_exception
-
+        void setDelegate(string&, function[void(vector[CxxEigenTriplet]&)], string&) except +translate_exception
 
 cdef extern from "cantera/cython/funcWrapper.h":
     # pyOverride is actually a templated function, but we have to specify the individual
@@ -77,6 +86,9 @@ cdef extern from "cantera/cython/funcWrapper.h":
     cdef function[int(string&, size_t)] pyOverride(PyObject*, int(PyFuncInfo&, string&, size_t))
     cdef function[int(size_t&, const string&)] pyOverride(
         PyObject*, int(PyFuncInfo&, size_t&, const string&))
+    cdef function[void(vector[CxxEigenTriplet]&)] pyOverride(
+        PyObject*, void(PyFuncInfo&, vector[CxxEigenTriplet]&))
+
 
 cdef extern from "cantera/base/ExtensionManager.h" namespace "Cantera":
     cdef cppclass CxxExtensionManager "Cantera::ExtensionManager":
diff --git a/interfaces/cython/cantera/delegator.pyx b/interfaces/cython/cantera/delegator.pyx
index c5b01bfd53..2166b245f7 100644
--- a/interfaces/cython/cantera/delegator.pyx
+++ b/interfaces/cython/cantera/delegator.pyx
@@ -9,7 +9,7 @@ from libc.stdlib cimport malloc
 from libc.string cimport strcpy
 
 from ._utils import CanteraError
-from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap
+from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap, get_triplet
 from .units cimport Units, UnitStack
 # from .reaction import ExtensibleRate, ExtensibleRateData
 from .reaction cimport (ExtensibleRate, ExtensibleRateData, CxxReaction,
@@ -244,6 +244,20 @@ cdef void callback_v_d_dp_dp(PyFuncInfo& funcInfo, size_array2 sizes, double arg
         funcInfo.setExceptionType(<PyObject*>exc_type)
         funcInfo.setExceptionValue(<PyObject*>exc_value)
 
+# Wrapper for void(vector<Eigen::Triplet>&)
+cdef void callback_v_vet(PyFuncInfo& funcInfo, vector[CxxEigenTriplet]& jac_vector) noexcept:
+    try:
+        python_trips = []
+        # convert vector to python object
+        (<object>funcInfo.func())(python_trips)
+        # add the triplets to the jacobian vector
+        for r, c, v in python_trips:
+            jac_vector.push_back(get_triplet(r, c, v))
+    except BaseException as e:
+        exc_type, exc_value = _sys.exc_info()[:2]
+        funcInfo.setExceptionType(<PyObject*>exc_type)
+        funcInfo.setExceptionValue(<PyObject*>exc_value)
+
 cdef int assign_delegates(obj, CxxDelegator* delegator) except -1:
     """
     Use methods defined in the Python class ``obj`` as delegates for the C++
@@ -366,6 +380,8 @@ cdef int assign_delegates(obj, CxxDelegator* delegator) except -1:
         elif callback == 'void(double,double*,double*)':
             delegator.setDelegate(cxx_name,
                 pyOverride(<PyObject*>method, callback_v_d_dp_dp), cxx_when)
+        elif callback == 'void(vector[CxxEigenTriplet]&)':
+            delegator.setDelegate(cxx_name, pyOverride(<PyObject*>method, callback_v_vet), cxx_when)
         else:
             raise ValueError("Don't know how to set delegates for functions "
                 f"with signature '{callback}'")
diff --git a/interfaces/cython/cantera/jacobians.pxd b/interfaces/cython/cantera/jacobians.pxd
index a5c6c45b18..3c0593b945 100644
--- a/interfaces/cython/cantera/jacobians.pxd
+++ b/interfaces/cython/cantera/jacobians.pxd
@@ -12,6 +12,8 @@ cdef extern from "cantera/numerics/SystemJacobian.h" namespace "Cantera":
         CxxSystemJacobian()
         string preconditionerSide()
         string type()
+        double gamma()
+        void setGamma(double) except +translate_exception
         void setPreconditionerSide(string) except +translate_exception
 
 cdef extern from "cantera/numerics/EigenSparseJacobian.h" namespace "Cantera":
@@ -44,20 +46,16 @@ cdef extern from "cantera/numerics/SystemJacobianFactory.h" namespace "Cantera":
 cdef class SystemJacobian:
     @staticmethod
     cdef wrap(shared_ptr[CxxSystemJacobian])
-    cdef set_cxx_object(self)
     cdef shared_ptr[CxxSystemJacobian] _base
 
 cdef class EigenSparseJacobian(SystemJacobian):
-    cdef set_cxx_object(self)
     cdef CxxEigenSparseJacobian* sparse_jac
 
 cdef class EigenSparseDirectJacobian(EigenSparseJacobian):
     pass
 
 cdef class AdaptivePreconditioner(EigenSparseJacobian):
-    cdef set_cxx_object(self)
     cdef CxxAdaptivePreconditioner* adaptive
 
 cdef class BandedJacobian(SystemJacobian):
-    cdef set_cxx_object(self)
     cdef CxxMultiJac* band_jac
diff --git a/interfaces/cython/cantera/jacobians.pyx b/interfaces/cython/cantera/jacobians.pyx
index bc951af06f..b6eaa3db74 100644
--- a/interfaces/cython/cantera/jacobians.pyx
+++ b/interfaces/cython/cantera/jacobians.pyx
@@ -48,7 +48,7 @@ cdef class SystemJacobian:
         jac.set_cxx_object()
         return jac
 
-    cdef set_cxx_object(self):
+    def set_cxx_object(self):
         pass
 
     property side:
@@ -63,6 +63,15 @@ cdef class SystemJacobian:
         def __set__(self, side):
             self._base.get().setPreconditionerSide(stringify(side))
 
+    property gamma:
+        """ Get/Set the value of gamma used in the expression P = (I - gamma * J).
+        """
+        def __get__(self):
+            return self._base.get().gamma()
+
+        def __set__(self, value):
+            self._base.get().setGamma(value)
+
 
 cdef class EigenSparseJacobian(SystemJacobian):
     """
@@ -72,7 +81,8 @@ cdef class EigenSparseJacobian(SystemJacobian):
 
     _type = "eigen-sparse"
 
-    cdef set_cxx_object(self):
+    def set_cxx_object(self):
+        super().set_cxx_object()
         self.sparse_jac = <CxxEigenSparseJacobian*>self._base.get()
 
     def print_contents(self):
@@ -108,7 +118,8 @@ cdef class AdaptivePreconditioner(EigenSparseJacobian):
     _type = "Adaptive"
     linear_solver_type = "GMRES"
 
-    cdef set_cxx_object(self):
+    def set_cxx_object(self):
+        super().set_cxx_object()
         self.adaptive = <CxxAdaptivePreconditioner*>self._base.get()
 
     property threshold:
@@ -179,5 +190,6 @@ cdef class BandedJacobian(SystemJacobian):
     _type = "banded-direct"
     linear_solver_type = "direct"
 
-    cdef set_cxx_object(self):
+    def set_cxx_object(self):
+        super().set_cxx_object()
         self.band_jac = <CxxMultiJac*>self._base.get()
diff --git a/interfaces/cython/cantera/reactor.pxd b/interfaces/cython/cantera/reactor.pxd
index ce460d10f8..9cf55612cb 100644
--- a/interfaces/cython/cantera/reactor.pxd
+++ b/interfaces/cython/cantera/reactor.pxd
@@ -5,6 +5,7 @@
 #distutils: language = c++
 
 from .ctcxx cimport *
+from .delegator cimport CxxEigenTriplet
 from .kinetics cimport *
 from .func1 cimport *
 from .jacobians cimport *
@@ -60,6 +61,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         void getState(double*) except +translate_exception
         CxxSparseMatrix jacobian() except +translate_exception
         CxxSparseMatrix finiteDifferenceJacobian() except +translate_exception
+        void buildJacobian(vector[CxxEigenTriplet]&) except +translate_exception
         void addSurface(CxxReactorSurface*) except +translate_exception
         void setAdvanceLimit(string&, double) except +translate_exception
         void addSensitivitySpeciesEnthalpy(size_t) except +translate_exception
@@ -210,6 +212,8 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
         void setPreconditioner(shared_ptr[CxxSystemJacobian] preconditioner)
         void setDerivativeSettings(CxxAnyMap&)
         CxxAnyMap solverStats() except +translate_exception
+        CxxSparseMatrix jacobian() except +translate_exception
+        CxxSparseMatrix finiteDifferenceJacobian() except +translate_exception
 
 cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera":
     cdef cppclass CxxReactorAccessor "Cantera::ReactorAccessor":
@@ -221,6 +225,7 @@ cdef extern from "cantera/zeroD/ReactorDelegator.h" namespace "Cantera":
         void setHeatRate(double)
         void restoreThermoState() except +translate_exception
         void restoreSurfaceState(size_t) except +translate_exception
+        void defaultEval(double time, double* LHS, double* RHS)
 
 
 ctypedef CxxReactorAccessor* CxxReactorAccessorPtr
diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
index b8715c98b1..dd2f381948 100644
--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -2,6 +2,8 @@
 # at https://cantera.org/license.txt for license and copyright information.
 
 import warnings
+import numpy as np
+from collections import defaultdict as _defaultdict
 import numbers as _numbers
 from cython.operator cimport dereference as deref
 import numpy as np
@@ -691,7 +693,8 @@ cdef class ExtensibleReactor(Reactor):
         'eval_surfaces': ('evalSurfaces', 'void(double*,double*,double*)'),
         'component_name': ('componentName', 'string(size_t)'),
         'component_index': ('componentIndex', 'size_t(string)'),
-        'species_index': ('speciesIndex', 'size_t(string)')
+        'species_index': ('speciesIndex', 'size_t(string)'),
+        'build_jacobian': ('buildJacobian', 'void(vector[CxxEigenTriplet]&)')
     }
 
     def __cinit__(self, *args, **kwargs):
@@ -750,6 +753,16 @@ cdef class ExtensibleReactor(Reactor):
         """
         self.accessor.restoreSurfaceState(n)
 
+    def default_eval(self, time, LHS, RHS):
+        """
+        Evaluation of the base reactors `eval` function to be used in `replace`
+        functions and maintain original functionality.
+        """
+        assert len(LHS) == self.n_vars and len(RHS) == self.n_vars
+        cdef np.ndarray[np.double_t, ndim=1, mode="c"] rhs = np.frombuffer(RHS)
+        cdef np.ndarray[np.double_t, ndim=1, mode="c"] lhs = np.frombuffer(LHS)
+        self.accessor.defaultEval(time, &lhs[0], &rhs[0])
+
 
 cdef class ExtensibleIdealGasReactor(ExtensibleReactor):
     """
@@ -2179,6 +2192,31 @@ cdef class ReactorNet:
         def __set__(self, settings):
             self.net.setDerivativeSettings(py_to_anymap(settings))
 
+    property jacobian:
+        """
+        Get the system Jacobian or an approximation thereof.
+
+        **Warning**: Depending on the particular implementation, this may return an
+        approximate Jacobian intended only for use in forming a preconditioner for
+        iterative solvers, excluding terms that would generate a fully-dense Jacobian.
+
+        **Warning**: This method is an experimental part of the Cantera API and may be
+        changed or removed without notice.
+        """
+        def __get__(self):
+            return get_from_sparse(self.net.jacobian(), self.n_vars, self.n_vars)
+
+    property finite_difference_jacobian:
+        """
+        Get the system Jacobian, calculated using a finite difference method.
+
+        **Warning:** this property is an experimental part of the Cantera API and
+        may be changed or removed without notice.
+        """
+        def __get__(self):
+            return get_from_sparse(self.net.finiteDifferenceJacobian(),
+                                   self.n_vars, self.n_vars)
+
     def draw(self, *, graph_attr=None, node_attr=None, edge_attr=None,
              heat_flow_attr=None, mass_flow_attr=None, moving_wall_edge_attr=None,
              surface_edge_attr=None, show_wall_velocity=True, print_state=False,
diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp
index ba0798eb21..7c49884823 100644
--- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp
+++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp
@@ -114,19 +114,21 @@ void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RH
     }
 }
 
-Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
+void IdealGasConstPressureMoleReactor::buildJacobian(
+    vector<Eigen::Triplet<double>>& jacVector)
 {
     if (m_nv == 0) {
         throw CanteraError("IdealGasConstPressureMoleReactor::jacobian",
                            "Reactor must be initialized first.");
     }
-    // clear former jacobian elements
-    m_jac_trips.clear();
     // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as
     // d (dot(omega)) / d c_j, it is later transformed appropriately.
     Eigen::SparseMatrix<double> dnk_dnj = m_kin->netProductionRates_ddCi();
     // species size that accounts for surface species
-    size_t ssize = m_nv - m_sidx;
+    size_t ssize = m_nsp;
+    for (auto surf : m_surfaces) {
+        ssize += surf->thermo()->nSpecies();
+    }
     // map derivatives from the surface chemistry jacobian
     // to the reactor jacobian
     if (!m_surfaces.empty()) {
@@ -169,7 +171,7 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
             if (static_cast<size_t>(it.row()) < m_nsp) {
                 it.valueRef() = it.value() + netProductionRates[it.row()] * molarVol;
             }
-            m_jac_trips.emplace_back(static_cast<int>(it.row() + m_sidx),
+            jacVector.emplace_back(static_cast<int>(it.row() + m_sidx),
                 static_cast<int>(it.col() + m_sidx), it.value());
         }
     }
@@ -198,9 +200,10 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
         for (size_t j = 0; j < m_nv; j++) {
             double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j];
             double ydotCurrent = rhsCurrent[j] / lhsCurrent[j];
-            m_jac_trips.emplace_back(static_cast<int>(j), 0,
+            jacVector.emplace_back(static_cast<int>(j), 0,
                                      (ydotPerturbed - ydotCurrent) / deltaTemp);
         }
+
         // d T_dot/dnj
         // allocate vectors for whole system
         Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize);
@@ -223,14 +226,16 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
         Eigen::VectorXd hk_dnkdnj_sums = dnk_dnj.transpose() * enthalpy;
         // Add derivatives to jac by spanning columns
         for (size_t j = 0; j < ssize; j++) {
-            m_jac_trips.emplace_back(0, static_cast<int>(j + m_sidx),
+            jacVector.emplace_back(0, static_cast<int>(j + m_sidx),
                 (specificHeat[j] * qdot - NCp * hk_dnkdnj_sums[j]) * denom);
         }
+
+        // build wall jacobian
+        buildWallJacobian(jacVector);
     }
-    // convert triplets to sparse matrix
-    Eigen::SparseMatrix<double> jac(m_nv, m_nv);
-    jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
-    return jac;
+
+    // build flow jacobian
+    buildFlowJacobian(jacVector);
 }
 
 size_t IdealGasConstPressureMoleReactor::componentIndex(const string& nm) const
@@ -287,4 +292,11 @@ double IdealGasConstPressureMoleReactor::lowerBound(size_t k) const {
     }
 }
 
+double IdealGasConstPressureMoleReactor::temperature_ddni(size_t index)
+{
+    // derivative of temperature transformed by ideal gas law
+    double n_total = m_mass / m_thermo->meanMolecularWeight();
+    return pressure() * m_vol / GasConstant / n_total;
+}
+
 }
diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp
index d7d39a7140..26f7026fb5 100644
--- a/src/zeroD/IdealGasMoleReactor.cpp
+++ b/src/zeroD/IdealGasMoleReactor.cpp
@@ -183,19 +183,20 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS)
     }
 }
 
-Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
+void IdealGasMoleReactor::buildJacobian(vector<Eigen::Triplet<double>>& jacVector)
 {
     if (m_nv == 0) {
         throw CanteraError("IdealGasMoleReactor::jacobian",
                            "Reactor must be initialized first.");
     }
-    // clear former jacobian elements
-    m_jac_trips.clear();
     // dnk_dnj represents d(dot(n_k)) / d (n_j) but is first assigned as
     // d (dot(omega)) / d c_j, it is later transformed appropriately.
     Eigen::SparseMatrix<double> dnk_dnj = m_kin->netProductionRates_ddCi();
     // species size that accounts for surface species
-    size_t ssize = m_nv - m_sidx;
+    size_t ssize = m_nsp;
+    for (auto surf : m_surfaces) {
+        ssize += surf->thermo()->nSpecies();
+    }
     // map derivatives from the surface chemistry jacobian
     // to the reactor jacobian
     if (!m_surfaces.empty()) {
@@ -215,7 +216,7 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
     // as it substantially reduces matrix sparsity
     for (int k = 0; k < dnk_dnj.outerSize(); k++) {
         for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
-            m_jac_trips.emplace_back(static_cast<int>(it.row() + m_sidx),
+            jacVector.emplace_back(static_cast<int>(it.row() + m_sidx),
                 static_cast<int>(it.col() + m_sidx), it.value());
         }
     }
@@ -243,9 +244,10 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
         for (size_t j = 0; j < m_nv; j++) {
             double ydotPerturbed = rhsPerturbed[j] / lhsPerturbed[j];
             double ydotCurrent = rhsCurrent[j] / lhsCurrent[j];
-            m_jac_trips.emplace_back(static_cast<int>(j), 0,
+            jacVector.emplace_back(static_cast<int>(j), 0,
                                      (ydotPerturbed - ydotCurrent) / deltaTemp);
         }
+
         // d T_dot/dnj
         Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize);
         Eigen::VectorXd internal_energy = Eigen::VectorXd::Zero(ssize);
@@ -272,14 +274,19 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
         Eigen::VectorXd uk_dnkdnj_sums = dnk_dnj.transpose() * internal_energy;
         // add derivatives to jacobian
         for (size_t j = 0; j < ssize; j++) {
-            m_jac_trips.emplace_back(0, static_cast<int>(j + m_sidx),
+            jacVector.emplace_back(0, static_cast<int>(j + m_sidx),
                 (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom);
         }
+        buildWallJacobian(jacVector);
     }
-    // convert triplets to sparse matrix
-    Eigen::SparseMatrix<double> jac(m_nv, m_nv);
-    jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
-    return jac;
+    buildFlowJacobian(jacVector);
+}
+
+double IdealGasMoleReactor::temperature_ddni(size_t index)
+{
+    // derivative of temperature transformed by ideal gas law
+    double n_total = m_mass / m_thermo->meanMolecularWeight();
+    return pressure() * m_vol / GasConstant / n_total;
 }
 
 }
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index bb6d8f5fae..88c00dcf19 100644
--- a/src/zeroD/Reactor.cpp
+++ b/src/zeroD/Reactor.cpp
@@ -55,6 +55,17 @@ void Reactor::setDerivativeSettings(AnyMap& settings)
     for (auto S : m_surfaces) {
         S->kinetics()->setDerivativeSettings(settings);
     }
+
+    // set reactor settings
+    bool force = settings.empty();
+    if (force || settings.hasKey("skip-walls")) {
+        m_jac_skip_walls = settings.getBool("skip-walls",
+            false);
+    }
+    if (force || settings.hasKey("skip-flow-devices")) {
+        m_jac_skip_flow_devices = settings.getBool("skip-flow-devices",
+            false);
+    }
 }
 
 void Reactor::setKinetics(Kinetics& kin)
@@ -662,4 +673,35 @@ void Reactor::setAdvanceLimit(const string& nm, const double limit)
     }
 }
 
+void Reactor::buildWallJacobian(vector<Eigen::Triplet<double>>& jacVector)
+{
+    if (!m_jac_skip_walls) {
+        for (size_t i = 0; i < m_wall.size(); i++) {
+            m_wall[i]->buildReactorJacobian(this, jacVector);
+        }
+    }
+}
+
+void Reactor::buildFlowJacobian(vector<Eigen::Triplet<double>>& jacVector)
+{
+    if (!m_jac_skip_flow_devices) {
+        for (size_t i = 0; i < m_outlet.size(); i++) {
+            m_outlet[i]->buildReactorJacobian(this, jacVector);
+        }
+
+        for (size_t i = 0; i <m_inlet.size(); i++) {
+            m_inlet[i]->buildReactorJacobian(this, jacVector);
+        }
+    }
+}
+
+Eigen::SparseMatrix<double> Reactor::jacobian() {
+        m_jac_trips.clear();
+        // Add before, during, after evals
+        buildJacobian(m_jac_trips);
+        // construct jacobian from vector
+        Eigen::SparseMatrix<double> jac(m_nv, m_nv);
+        jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end());
+        return jac;
+    }
 }
diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp
index d0e73adc74..dbc22d9417 100644
--- a/src/zeroD/ReactorNet.cpp
+++ b/src/zeroD/ReactorNet.cpp
@@ -160,7 +160,26 @@ void ReactorNet::initialize()
                 " a single ReactorNet will be an error after Cantera 3.2.", shared);
         }
     }
-
+    // Create walls and flow devices sets
+    for (auto r : m_reactors) {
+        // walls
+        if (!m_jac_skip_walls) {
+            for (size_t i = 0; i < r->nWalls(); i++) {
+                m_walls.insert(&(r->wall(i)));
+            }
+        }
+        // flow devices
+        if (!m_jac_skip_flow_devices) {
+            // outlets
+            for (size_t i = 0; i < r->nOutlets(); i++) {
+                m_flow_devices.insert(&(r->outlet(i)));
+            }
+            // inlets
+            for (size_t i = 0; i < r->nInlets(); i++) {
+                m_flow_devices.insert(&(r->inlet(i)));
+            }
+        }
+    }
     m_ydot.resize(m_nv,0.0);
     m_yest.resize(m_nv,0.0);
     m_advancelimits.resize(m_nv,-1.0);
@@ -691,12 +710,32 @@ size_t ReactorNet::registerSensitivityParameter(
     return m_sens_params.size() - 1;
 }
 
+size_t ReactorNet::globalStartIndex(ReactorBase* curr_reactor) {
+        for (size_t i = 0; i < m_reactors.size(); i++) {
+            if (curr_reactor == m_reactors[i]) {
+                return m_start[i];
+            }
+        }
+        throw CanteraError("ReactorNet::globalStartIndex: ",
+                curr_reactor->name(), " not found in network.");
+    }
+
 void ReactorNet::setDerivativeSettings(AnyMap& settings)
 {
     // Apply given settings to all reactors
     for (size_t i = 0; i < m_reactors.size(); i++) {
         m_reactors[i]->setDerivativeSettings(settings);
     }
+    // set network settings
+    bool force = settings.empty();
+    if (force || settings.hasKey("skip-walls")) {
+        m_jac_skip_walls = settings.getBool("skip-walls",
+            false);
+    }
+    if (force || settings.hasKey("skip-flow-devices")) {
+        m_jac_skip_flow_devices = settings.getBool("skip-flow-devices",
+            false);
+    }
 }
 
 AnyMap ReactorNet::solverStats() const
@@ -740,19 +779,16 @@ void ReactorNet::preconditionerSetup(double t, double* y, double gamma)
     vector<double> yCopy(m_nv);
     // Get state of reactor
     getState(yCopy.data());
-    // transform state based on preconditioner rules
+    // Transform state based on preconditioner rules
     precon->stateAdjustment(yCopy);
-    // update network with adjusted state
+    // Update network with adjusted state
     updateState(yCopy.data());
-    // Get jacobians and give elements to preconditioners
-    for (size_t i = 0; i < m_reactors.size(); i++) {
-        Eigen::SparseMatrix<double> rJac = m_reactors[i]->jacobian();
-        for (int k=0; k<rJac.outerSize(); ++k) {
-            for (Eigen::SparseMatrix<double>::InnerIterator it(rJac, k); it; ++it) {
-                precon->setValue(it.row() + m_start[i], it.col() + m_start[i],
-                    it.value());
-            }
-        }
+    // Create jacobian triplet vector
+    vector<Eigen::Triplet<double>> jacVector;
+    buildJacobian(jacVector);
+    // Add to preconditioner with offset
+    for (auto it : jacVector) {
+        precon->setValue(it.row(), it.col(), it.value());
     }
     // post reactor setup operations
     precon->updatePreconditioner();
@@ -781,6 +817,87 @@ void ReactorNet::checkPreconditionerSupported() const {
     }
 }
 
+void ReactorNet::buildJacobian(vector<Eigen::Triplet<double>>& jacVector)
+{
+    // network must be initialized for the jacobian
+    if (!m_init) {
+        initialize();
+    }
+    // Create jacobian triplet vector
+    vector<size_t> jstarts;
+    // Get jacobians and give elements to preconditioners
+    jstarts.push_back(jacVector.size());
+    for (size_t i = 0; i < m_reactors.size(); i++) {
+        m_reactors[i]->buildJacobian(jacVector);
+        jstarts.push_back(jacVector.size());
+    }
+    // Add to preconditioner with offset
+    for (size_t i=0; i < m_reactors.size(); i++) {
+        for (size_t j = jstarts[i]; j < jstarts[i+1]; j++) {
+            auto it = jacVector[j];
+            auto newTrip = Eigen::Triplet<double>(it.row() + m_start[i], it.col()
+                + m_start[i], it.value());
+            jacVector[j] = newTrip;
+        }
+    }
+
+    // loop through all connections and then set them found so calculations are not
+    // repeated
+    for (auto r : m_reactors) {
+        // walls
+        for (const auto& wall : m_walls) {
+            wall->buildNetworkJacobian(jacVector);
+        }
+        // flow devices
+        for (const auto& flow_device : m_flow_devices) {
+            flow_device->buildNetworkJacobian(jacVector);
+        }
+    }
+}
+
+Eigen::SparseMatrix<double> ReactorNet::finiteDifferenceJacobian()
+{
+    // network must be initialized for the jacobian
+    if (! m_init) {
+        initialize();
+    }
+
+    // allocate jacobian triplet vector
+    vector<Eigen::Triplet<double>> jac_trips;
+
+    // Get the current state
+    Eigen::ArrayXd yCurrent(m_nv);
+    getState(yCurrent.data());
+
+    Eigen::ArrayXd yPerturbed = yCurrent;
+    Eigen::ArrayXd ydotCurrent(m_nv), ydotPerturbed(m_nv);
+
+    eval(m_time, yCurrent.data(), ydotCurrent.data(), m_sens_params.data());
+    double rel_perturb = std::sqrt(std::numeric_limits<double>::epsilon());
+
+    for (size_t j = 0; j < m_nv; j++) {
+        yPerturbed = yCurrent;
+        double delta_y = std::max(std::abs(yCurrent[j]), 1000 * m_atols) * rel_perturb;
+        yPerturbed[j] += delta_y;
+        ydotPerturbed = 0;
+        eval(m_time, yPerturbed.data(), ydotPerturbed.data(), m_sens_params.data());
+        // d ydot_i/dy_j
+        for (size_t i = 0; i < m_nv; i++) {
+            if (ydotCurrent[i] != ydotPerturbed[i]) {
+                jac_trips.emplace_back(
+                    static_cast<int>(i), static_cast<int>(j),
+                    (ydotPerturbed[i] - ydotCurrent[i]) / delta_y);
+            }
+        }
+    }
+    updateState(yCurrent.data());
+
+    Eigen::SparseMatrix<double> jac(m_nv, m_nv);
+    jac.setFromTriplets(jac_trips.begin(), jac_trips.end());
+    return jac;
+}
+
+
 SteadyReactorSolver::SteadyReactorSolver(ReactorNet* net, double* x0)
     : m_net(net)
 {
diff --git a/src/zeroD/Wall.cpp b/src/zeroD/Wall.cpp
index 9d638c3281..57a6bdec3a 100644
--- a/src/zeroD/Wall.cpp
+++ b/src/zeroD/Wall.cpp
@@ -6,6 +6,8 @@
 #include "cantera/base/stringUtils.h"
 #include "cantera/numerics/Func1.h"
 #include "cantera/zeroD/Wall.h"
+#include "cantera/thermo/ThermoPhase.h"
+#include "cantera/zeroD/ReactorNet.h"
 
 namespace Cantera
 {
@@ -93,4 +95,66 @@ double Wall::heatRate()
     return q1;
 }
 
+
+void Wall::buildReactorJacobian(ReactorBase* r,
+    vector<Eigen::Triplet<double>>& jacVector)
+{
+    // get derivative of heat transfer for both reactors
+    vector<Eigen::Triplet<double>> network;
+    size_t nsp = r->contents().nSpecies();
+    size_t sidx = r->speciesOffset();
+    size_t eidx = r->energyIndex();
+    // define a scalar for direction based on left and right
+    double direction = (r == m_left) ? 1.0 : -1.0;
+    // elements within the current reactor
+    // find dQdni for the current reactor w.r.t current reactor
+    for (size_t i = sidx; i < nsp + sidx; i++) {
+        double dQdni = m_rrth * m_area * direction * r->temperature_ddni(i);
+        dQdni += m_emiss * m_area * direction * r->temperature_ddni(i) * 4
+                 * pow(r->temperature(), 3);
+        jacVector.emplace_back(eidx, i, dQdni);
+    }
+}
+
+void Wall::buildNetworkJacobian(vector<Eigen::Triplet<double>>& jacVector)
+{
+    // No interdependent terms for reservoirs
+    if (m_right->type() == "Reservoir" || m_left->type() == "Reservoir") {
+        return;
+    }
+    // get derivatives for inter-dependent reactor terms
+    //variables for the right side
+    vector<Eigen::Triplet<double>> network;
+    size_t r_nsp = m_right->contents().nSpecies();
+    size_t r_sidx = m_right->speciesOffset();
+    size_t r_net = m_right->network().globalStartIndex(m_right);
+    size_t r_eidx = m_right->energyIndex();
+
+    // variables for the left side
+    size_t l_nsp = m_left->contents().nSpecies();
+    size_t l_sidx = m_left->speciesOffset();
+    size_t l_net = m_left->network().globalStartIndex(m_left);
+    size_t l_eidx = m_left->energyIndex();
+
+    if (m_right->energyEnabled()) {
+        // find dQdni for the right reactor w.r.t left reactor
+        for (size_t i = l_sidx; i < l_sidx + l_nsp; i++) {
+            double dQdni = m_rrth * m_area * m_left->temperature_ddni(i);
+            dQdni += m_emiss * m_area * m_left->temperature_ddni(i) * 4
+                 * pow(m_left->temperature(), 3);
+            jacVector.emplace_back(r_eidx + r_net, i + l_net, dQdni);
+        }
+    }
+
+    if (m_left->energyEnabled()) {
+        // find dQdni for the left reactor w.r.t right reactor
+        for (size_t i = r_sidx; i < r_sidx + r_nsp; i++) {
+            double dQdni = - m_rrth * m_area * m_right->temperature_ddni(i);
+            dQdni -= m_emiss * m_area * m_right->temperature_ddni(i) * 4
+                 * pow(m_right->temperature(), 3);
+            jacVector.emplace_back(l_eidx + l_net, i + r_net, dQdni);
+        }
+    }
+}
+
 }
diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py
index 33c882194f..cb4ceebdb2 100644
--- a/test/python/test_reactor.py
+++ b/test/python/test_reactor.py
@@ -185,6 +185,61 @@ def test_finite_difference_jacobian(self):
             if name in constant:
                 assert all(J[i, species_start:] == 0), (i, name)
 
+    def test_network_finite_difference_jacobian(self):
+        self.make_reactors(T1=900, P1=101325, X1="H2:0.4, O2:0.4, N2:0.2")
+        k1H2 = self.gas1.species_index("H2")
+        k2H2 = self.gas1.species_index("H2")
+        while self.r1.thermo.X[k1H2] > 0.3 or self.r2.thermo.X[k2H2] > 0.3:
+            self.net.step()
+
+        J = self.net.finite_difference_jacobian
+        assert J.shape == (self.net.n_vars, self.net.n_vars)
+
+        # state variables that should be constant, depending on reactor type
+        constant = {"mass", "volume", "int_energy", "enthalpy", "pressure"}
+        variable = {"temperature"}
+        for i in range(3):
+            name = self.r1.component_name(i)
+            if name in constant:
+                assert all(J[i,:] == 0), (i, name)
+            elif name in variable:
+                assert any(J[i,:] != 0)
+            # check in second reactor
+            name = self.r2.component_name(i)
+            if name in constant:
+                assert all(J[i + self.r1.n_vars,:] == 0), (i, name)
+            elif name in variable:
+                assert any(J[i + self.r1.n_vars,:] != 0)
+
+        # Disabling energy equation should zero these terms
+        self.r1.energy_enabled = False
+        self.r2.energy_enabled = False
+        J = self.net.finite_difference_jacobian
+        for i in range(3):
+            name = self.r1.component_name(i)
+            if name == "temperature":
+                assert all(J[i,:] == 0)
+            name = self.r2.component_name(i)
+            if name == "temperature":
+                assert all(J[i + self.r1.n_vars,:] == 0)
+
+        # Disabling species equations should zero these terms
+        self.r1.energy_enabled = True
+        self.r1.chemistry_enabled = False
+        self.r2.energy_enabled = True
+        self.r2.chemistry_enabled = False
+        J = self.net.finite_difference_jacobian
+        constant = set(self.gas1.species_names + self.gas2.species_names)
+        r1_species_start = self.r1.component_index(self.gas1.species_name(0))
+        r2_species_start = self.r2.component_index(self.gas2.species_name(0))
+        for i in range(self.r1.n_vars):
+            name = self.r1.component_name(i)
+            if name in constant:
+                assert all(J[i, r1_species_start:] == 0), (i, name)
+            name = self.r2.component_name(i)
+            if name in constant:
+                assert all(J[i + self.r1.n_vars, (r2_species_start + self.r1.n_vars):] == 0), (i, name)
+
     def test_timestepping(self):
         self.make_reactors()
 
@@ -1462,7 +1517,7 @@ def create_reactors(self, **kwargs):
         self.precon = ct.AdaptivePreconditioner()
         self.net2.preconditioner = self.precon
         self.net2.derivative_settings = {"skip-third-bodies":True, "skip-falloff":True,
-            "skip-coverage-dependence":True}
+            "skip-coverage-dependence":True, "skip-flow-devices": True}
 
     def test_get_solver_type(self):
         self.create_reactors()
@@ -1470,11 +1525,70 @@ def test_get_solver_type(self):
         self.net2.initialize()
         assert self.net2.linear_solver_type == "GMRES"
 
+    def test_mass_flow_jacobian(self):
+        self.create_reactors(add_mdot=True)
+        # reset derivative settings
+        self.net2.derivative_settings = {"skip-third-bodies":True, "skip-falloff":True,
+            "skip-coverage-dependence":True, "skip-flow-devices": False}
+
+        with pytest.raises(NotImplementedError, match="MassFlowController::buildReactorJacobian"):
+            J = self.net2.jacobian
+
+        with pytest.raises(NotImplementedError, match="MassFlowController::buildReactorJacobian"):
+            J = self.r2.jacobian
+
+    @pytest.mark.xfail
+    def test_heat_transfer_network(self):
+        # create first reactor
+        gas1 = ct.Solution("h2o2.yaml", "ohmech")
+        gas1.TPX = 600, ct.one_atm, "O2:1.0"
+        r1 = self.reactorClass(gas1)
+
+        # create second reactor
+        gas2 = ct.Solution("h2o2.yaml", "ohmech")
+        gas2.TPX = 300, ct.one_atm, "O2:1.0"
+        r2 = ct.reactorClass(gas2)
+
+        # create wall
+        U = 2.
+        A = 3.0
+        w = ct.Wall(r1, r2, U=U, A=A)
+        net = ct.ReactorNet([r1,])
+        jac = net.jacobian
+        fd_jac = net.finite_difference_jacobian
+        for i in range(jac.shape[0]):
+            for j in range(jac.shape[1]):
+                assert np.isclose(jac[i, j], fd_jac[i, j])
 
 class TestIdealGasMoleReactor(TestMoleReactor):
     reactorClass = ct.IdealGasMoleReactor
     test_preconditioner_unsupported = None
 
+    @pytest.mark.xfail
+    def test_heat_transfer_network(self):
+        # create first reactor
+        gas1 = ct.Solution("h2o2.yaml", "ohmech")
+        gas1.TPX = 600, ct.one_atm, "O2:1.0"
+        r1 = self.reactorClass(gas1)
+
+        # create second reactor
+        gas2 = ct.Solution("h2o2.yaml", "ohmech")
+        gas2.TPX = 300, ct.one_atm, "O2:1.0"
+        r2 = self.reactorClass(gas2)
+
+        # create wall
+        U = 2.0
+        A = 3.0
+        w = ct.Wall(r1, r2, U=U, A=A)
+        net = ct.ReactorNet([r1, r2])
+        jac = net.jacobian
+        fd_jac = net.finite_difference_jacobian
+        # check for values
+        for i in range(jac.shape[0]):
+            for j in range(jac.shape[1]):
+                assert np.isclose(jac[i, j], fd_jac[i, j])
+
+
     def test_adaptive_precon_integration(self):
         # Network one with non-mole reactor
         T0 = 900
@@ -3269,6 +3383,109 @@ def deltaC():
         assert r1.phase["H2"].Y[0] == approx(0.13765976, rel=1e-6)
         assert r2.phase["O2"].Y[0] == approx(0.94617029, rel=1e-6)
 
+    def test_after_jacobian(self):
+        class AfterJacobianReactor(ct.ExtensibleIdealGasMoleReactor):
+            def __init__(self, *args, neighbor, **kwargs):
+                super().__init__(*args, **kwargs)
+                self.v_wall = 0
+                self.k_wall = 1e-5
+                self.neighbor = neighbor
+
+            def after_initialize(self, t0):
+                self.n_vars += 1
+                self.i_wall = self.n_vars - 1
+
+            def after_get_state(self, y):
+                y[self.i_wall] = self.v_wall
+
+            def after_update_state(self, y):
+                self.v_wall = y[self.i_wall]
+                self.walls[0].velocity = self.v_wall
+
+            def after_eval(self, t, LHS, RHS):
+                # Extra equation is d(v_wall)/dt = k * delta P
+                a = self.k_wall * (self.thermo.P - self.neighbor.thermo.P)
+                RHS[self.i_wall] = a
+
+            def before_component_index(self, name):
+                if name == 'v_wall':
+                    return self.i_wall
+
+            def before_component_name(self, i):
+                if i == self.i_wall:
+                    return 'v_wall'
+
+            def after_build_jacobian(self, jac_vector):
+                jac_vector.append((self.i_wall, self.i_wall, 1e20))
+
+        self.gas.TP = 300, ct.one_atm
+        res = ct.Reservoir(self.gas)
+        self.gas.TP = 300, 2 * ct.one_atm
+        r = AfterJacobianReactor(self.gas, neighbor=res)
+        w = ct.Wall(r, res)
+        net = ct.ReactorNet([r])
+        precon = ct.AdaptivePreconditioner()
+        net.preconditioner = precon
+        net.step()
+        # test that jacobian wall element is hard coded value
+        jac = r.jacobian
+        assert jac[r.i_wall, r.i_wall] == 1e20
+        pmat = precon.matrix
+        assert pmat[r.i_wall, r.i_wall] == (1 - precon.gamma * 1e20)
+
+    def test_before_jacobian(self):
+        class BeforeJacobianReactor(ct.ExtensibleIdealGasMoleReactor):
+
+            def before_build_jacobian(self, jac_vector):
+                jac_vector.append((0, 0, 1e10))
+
+        self.gas.TP = 300, ct.one_atm
+        r = BeforeJacobianReactor(self.gas)
+        net = ct.ReactorNet([r])
+        net.preconditioner = ct.AdaptivePreconditioner()
+        net.step()
+        # test that jacobian wall element is hard coded value
+        jac = r.jacobian
+        assert jac[0, 0] == 1e10
+
+    def test_replace_jacobian(self):
+        class ReplaceJacobianReactor(ct.ExtensibleIdealGasMoleReactor):
+
+            def replace_build_jacobian(self, jac_vector):
+                jac_vector.append((0, 0, 0))
+
+        self.gas.TP = 300, ct.one_atm
+        r = ReplaceJacobianReactor(self.gas)
+        net = ct.ReactorNet([r])
+        net.preconditioner = ct.AdaptivePreconditioner()
+        net.step()
+        # test that jacobian wall element is hard coded value
+        jac = r.jacobian
+        assert np.sum(jac) == 0
+
+    def test_replace_with_default_eval(self):
+        class ReplaceEvalReactor(ct.ExtensibleIdealGasConstPressureMoleReactor):
+
+            def replace_eval(self, t, LHS, RHS):
+                self.default_eval(t, LHS, RHS)
+
+        # setup thermo object
+        gas = ct.Solution("h2o2.yaml", "ohmech")
+        gas.set_equivalence_ratio(0.5, "H2:1.0", "O2:1.0")
+        gas.equilibrate("HP")
+        # replacement reactor
+        r = ReplaceEvalReactor(gas)
+        r.volume = 1.0
+        # default reactor
+        rstd = ct.IdealGasConstPressureMoleReactor(gas)
+        rstd.volume = r.volume
+        # network of both reactors
+        net = ct.ReactorNet([r, rstd])
+        net.preconditioner = ct.AdaptivePreconditioner()
+        net.advance_to_steady_state()
+        # reactors should have the same solution because the default is used
+        assert r.get_state() == approx(rstd.get_state())
+
 
 class TestSteadySolver:
     @pytest.mark.parametrize("reactor_class",
diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp
index be07f750a6..46ef4b9162 100644
--- a/test/zeroD/test_zeroD.cpp
+++ b/test/zeroD/test_zeroD.cpp
@@ -307,6 +307,106 @@ TEST(AdaptivePreconditionerTests, test_precon_solver_stats)
     EXPECT_GE(stats["nonlinear_conv_fails"].asInt(), 0);
 }
 
+TEST(JacobianTests, test_wall_jacobian_build)
+{
+    // create first reactor
+    auto sol1 = newSolution("h2o2.yaml");
+    sol1->thermo()->setState_TPY(1000.0, OneAtm, " O2:1.0");
+    auto reactor1 = make_shared<IdealGasMoleReactor>(sol1);
+    reactor1->setInitialVolume(1.0);
+    // create second reactor
+    auto sol2 = newSolution("h2o2.yaml");
+    sol2->thermo()->setState_TPY(900.0, OneAtm, " O2:1.0");
+    auto reactor2 = make_shared<IdealGasConstPressureMoleReactor>(sol2);
+    reactor2->setInitialVolume(1.0);
+    // create the wall
+    Wall w(reactor1, reactor2);
+    w.setArea(2.0);
+    w.setHeatTransferCoeff(3.0);
+    // setup reactor network and integrate
+    vector<shared_ptr<ReactorBase>> reactors = {reactor1, reactor2};
+    ReactorNet network(reactors);
+    network.initialize();
+    // create jacobian the size of network
+    Eigen::SparseMatrix<double> wallJacMat;
+    wallJacMat.resize(network.neq(), network.neq());
+    // manually get wall jacobian elements
+    vector<Eigen::Triplet<double>>  wallJac;
+    // build jac for reactor 1 wall only
+    w.buildReactorJacobian(reactor1.get(), wallJac);
+    wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end());
+    // check that wall jacobian forms correct value
+    double v1 = sol1->thermo()->temperature() * w.area() * w.getHeatTransferCoeff();
+    for (int k = 0; k < wallJacMat.outerSize(); k++) {
+        for (Eigen::SparseMatrix<double>::InnerIterator it(wallJacMat, k); it; ++it) {
+            EXPECT_DOUBLE_EQ(it.value(), v1);
+            EXPECT_EQ(it.row(), 0); // check that it is the first row
+            EXPECT_GE(it.col(), reactor1->speciesOffset());
+            EXPECT_LT(it.col(), reactor1->neq());
+        }
+    }
+    // build jac for reactor 2 wall only
+    wallJac.clear();
+    w.buildReactorJacobian(reactor2.get(), wallJac);
+    wallJacMat.setZero();
+    wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end());
+    // check that wall jacobian forms correct value
+    double v2 = sol2->thermo()->temperature() * w.area() * w.getHeatTransferCoeff();
+    for (int k = 0; k < wallJacMat.outerSize(); k++) {
+        for (Eigen::SparseMatrix<double>::InnerIterator it(wallJacMat, k); it; ++it) {
+            EXPECT_DOUBLE_EQ(it.value(), -v2);
+            EXPECT_EQ(it.row(), 0); // check that it is the first row
+            EXPECT_GE(it.col(), reactor2->speciesOffset());
+            EXPECT_LT(it.col(), reactor2->neq());
+        }
+    }
+    // build jac for network terms
+    wallJac.clear();
+    w.buildNetworkJacobian(wallJac);
+    wallJacMat.setZero();
+    wallJacMat.setFromTriplets(wallJac.begin(), wallJac.end());
+    // check appropriate values
+    // double tol = 1e-8;
+    for (int k = 0; k < wallJacMat.outerSize(); k++) {
+        for (Eigen::SparseMatrix<double>::InnerIterator it(wallJacMat, k); it; ++it) {
+            if (it.value() < 0) {
+                EXPECT_DOUBLE_EQ(it.value(), -5400.0);
+                EXPECT_EQ(it.row(), 0); // check that it is the first row
+                EXPECT_GE(it.col(), reactor1->neq() + reactor2->speciesOffset());
+                EXPECT_LT(it.col(), reactor1->neq() + reactor2->neq());
+            } else {
+                EXPECT_DOUBLE_EQ(it.value(), 6000.0);
+                EXPECT_EQ(it.row(), reactor1->neq()); // check that it is the first row
+                EXPECT_GE(it.col(), reactor2->speciesOffset());
+                EXPECT_LT(it.col(), reactor1->neq());
+            }
+        }
+    }
+}
+
+TEST(JacobianTests, test_flow_jacobian_not_implemented)
+{
+    // create reservoir reactor
+    auto sol = newSolution("h2o2.yaml");
+    sol->thermo()->setState_TPY(1000.0, OneAtm, "O2:1.0");
+    auto res = make_shared<Reservoir>(sol);
+    // create reactor
+    auto reactor = make_shared<IdealGasConstPressureMoleReactor>(sol);
+    reactor->setInitialVolume(1.0);
+    // create the flow device
+    MassFlowController mfc(res, reactor);
+    mfc.setMassFlowCoeff(1.0);
+    // setup reactor network and integrate
+    ReactorNet network(reactor);
+    network.initialize();
+    // manually get wall jacobian elements
+    vector<Eigen::Triplet<double>>  flowJac;
+    // expect errors from building jacobians
+    EXPECT_THROW(mfc.buildReactorJacobian(reactor.get(), flowJac), NotImplementedError);
+    // check the jacobian calculated flag and throw/catch errors accordingly
+    EXPECT_THROW(mfc.buildNetworkJacobian(flowJac), NotImplementedError);
+}
+
 int main(int argc, char** argv)
 {
     printf("Running main() from test_zeroD.cpp\n");