diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
index de221c6c5b..9e5ef75ec7 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
@@ -153,7 +153,7 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omegaI
 template<class ReactionThermo, class ThermoType>
 void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
 (
-    const scalar time,
+    const scalar t,
     const scalarField& c,
     const label li,
     scalarField& dcdt
@@ -186,15 +186,15 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
     }
     cp /= rho;
 
-    scalar dT = 0;
+    scalar dTdt = 0;
     for (label i = 0; i < nSpecie_; i++)
     {
         const scalar hi = specieThermos_[i].ha(p, T);
-        dT += hi*dcdt[i];
+        dTdt += hi*dcdt[i];
     }
-    dT /= rho*cp;
+    dTdt /= rho*cp;
 
-    dcdt[nSpecie_] = -dT;
+    dcdt[nSpecie_] = -dTdt;
 
     // dp/dt = ...
     dcdt[nSpecie_ + 1] = 0;
@@ -257,6 +257,11 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
     }
     dTdt /= -cpMean; // K/s
 
+    dcdt[nSpecie_] = dTdt;
+
+    // dp/dt = ...
+    dcdt[nSpecie_ + 1] = 0;
+
     for (label i = 0; i < nSpecie_; i++)
     {
         J(nSpecie_, i) = 0;
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
index 0c9d8214f9..37f9ef36f1 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
@@ -71,9 +71,6 @@ class StandardChemistryModel
 
 protected:
 
-    typedef ThermoType thermoType;
-
-
     // Protected data
 
         //- Reference to the field of specie mass fractions
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
index 6956974997..ba49ff69a9 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
@@ -612,8 +612,8 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::solve
             c0[i] = c[i];
             phiq[i] = this->Y()[i][celli];
         }
-        phiq[this->nSpecie()]=Ti;
-        phiq[this->nSpecie() + 1]=pi;
+        phiq[this->nSpecie()] = Ti;
+        phiq[this->nSpecie() + 1] = pi;
         if (tabulation_->variableTimeStep())
         {
             phiq[this->nSpecie() + 2] = deltaT[celli];
diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
index aab1ae27c8..32ab244a08 100644
--- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
+++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
@@ -25,8 +25,6 @@ License
 
 #include "EulerImplicit.H"
 #include "addToRunTimeSelectionTable.H"
-#include "simpleMatrix.H"
-#include "Reaction.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -39,8 +37,10 @@ Foam::EulerImplicit<ChemistryModel>::EulerImplicit
     chemistrySolver<ChemistryModel>(thermo),
     coeffsDict_(this->subDict("EulerImplicitCoeffs")),
     cTauChem_(coeffsDict_.lookup<scalar>("cTauChem")),
-    eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")),
-    cTp_(this->nEqns())
+    cTp_(this->nEqns()),
+    R_(this->nEqns()),
+    J_(this->nEqns()),
+    E_(this->nEqns() - 2)
 {}
 
 
@@ -53,41 +53,6 @@ Foam::EulerImplicit<ChemistryModel>::~EulerImplicit()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ChemistryModel>
-void Foam::EulerImplicit<ChemistryModel>::updateRRInReactionI
-(
-    const label index,
-    const scalar pr,
-    const scalar pf,
-    const scalar corr,
-    const label lRef,
-    const label rRef,
-    const scalar p,
-    const scalar T,
-    simpleMatrix<scalar>& RR
-) const
-{
-    const Reaction<typename ChemistryModel::thermoType>& R =
-        this->reactions_[index];
-
-    forAll(R.lhs(), s)
-    {
-        const label si = R.lhs()[s].index;
-        const scalar sl = R.lhs()[s].stoichCoeff;
-        RR[si][rRef] -= sl*pr*corr;
-        RR[si][lRef] += sl*pf*corr;
-    }
-
-    forAll(R.rhs(), s)
-    {
-        const label si = R.rhs()[s].index;
-        const scalar sr = R.rhs()[s].stoichCoeff;
-        RR[si][lRef] -= sr*pf*corr;
-        RR[si][rRef] += sr*pr*corr;
-    }
-}
-
-
 template<class ChemistryModel>
 void Foam::EulerImplicit<ChemistryModel>::solve
 (
@@ -100,101 +65,71 @@ void Foam::EulerImplicit<ChemistryModel>::solve
 ) const
 {
     const label nSpecie = this->nSpecie();
-    simpleMatrix<scalar> RR(nSpecie, 0, 0);
 
-    for (label i=0; i<nSpecie; i++)
+    // Map the composition, temperature and pressure into cTp
+    for (int i=0; i<nSpecie; i++)
     {
-        c[i] = max(0, c[i]);
+        cTp_[i] = max(0, c[i]);
     }
+    cTp_[nSpecie] = T;
+    cTp_[nSpecie + 1] = p;
 
-    // Calculate the absolute enthalpy
-    const scalar cTot = sum(c);
-    typename ChemistryModel::thermoType mixture
-    (
-        (this->specieThermos_[0].W()*c[0])*this->specieThermos_[0]
-    );
-    for (label i=1; i<nSpecie; i++)
-    {
-        mixture += (this->specieThermos_[i].W()*c[i])*this->specieThermos_[i];
-    }
-    const scalar ha = mixture.Ha(p, T);
-    const scalar deltaTEst = min(deltaT, subDeltaT);
-
-    forAll(this->reactions(), i)
-    {
-        scalar pf, cf, pr, cr;
-        label lRef, rRef;
-
-        const scalar omegai
-        (
-            this->omegaI(i, p, T, c, li, pf, cf, lRef, pr, cr, rRef)
-        );
-
-        scalar corr = 1;
-        if (eqRateLimiter_)
-        {
-            if (omegai < 0)
-            {
-                corr = 1/(1 + pr*deltaTEst);
-            }
-            else
-            {
-                corr = 1/(1 + pf*deltaTEst);
-            }
-        }
-
-        updateRRInReactionI(i, pr, pf, corr, lRef, rRef, p, T, RR);
-    }
+    // Calculate the reaction rate and Jacobian
+    this->jacobian(0, cTp_, li, R_, J_);
 
     // Calculate the stable/accurate time-step
     scalar tMin = great;
+    const scalar cTot = sum(c);
 
     for (label i=0; i<nSpecie; i++)
     {
-        scalar d = 0;
-        for (label j=0; j<nSpecie; j++)
+        if (R_[i] < -small)
         {
-            d -= RR(i, j)*c[j];
-        }
-
-        if (d < -small)
-        {
-            tMin = min(tMin, -(c[i] + small)/d);
+            tMin = min(tMin, -(cTp_[i] + small)/R_[i]);
         }
         else
         {
-            d = max(d, small);
-            const scalar cm = max(cTot - c[i], 1e-5);
-            tMin = min(tMin, cm/d);
+            tMin = min
+            (
+                tMin,
+                max(cTot - cTp_[i], 1e-5)/max(R_[i], small)
+            );
         }
     }
 
     subDeltaT = cTauChem_*tMin;
     deltaT = min(deltaT, subDeltaT);
 
-    // Add the diagonal and source contributions from the time-derivative
+    // Assemble the Euler implicit matrix for the composition
+    scalarField& source = E_.source();
     for (label i=0; i<nSpecie; i++)
     {
-        RR(i, i) += 1/deltaT;
-        RR.source()[i] = c[i]/deltaT;
+        E_(i, i) = 1/deltaT - J_(i, i);
+        source[i] = R_[i] + E_(i, i)*cTp_[i];
+
+        for (label j=0; j<nSpecie; j++)
+        {
+            if (i != j)
+            {
+                E_(i, j) = -J_(i, j);
+                source[i] += E_(i, j)*cTp_[j];
+            }
+        }
     }
 
     // Solve for the new composition
-    c = RR.LUsolve();
+    scalarField::subField(cTp_, nSpecie) = E_.LUsolve();
 
-    // Limit the composition
+    // Limit the composition and transfer back into c
     for (label i=0; i<nSpecie; i++)
     {
-        c[i] = max(0, c[i]);
+        c[i] = max(0, cTp_[i]);
     }
 
-    // Update the temperature
-    mixture = (this->specieThermos_[0].W()*c[0])*this->specieThermos_[0];
-    for (label i=1; i<nSpecie; i++)
-    {
-        mixture += (this->specieThermos_[i].W()*c[i])*this->specieThermos_[i];
-    }
-    T = mixture.THa(ha, p, T);
+    // Euler explicit integrate the temperature.
+    // Separating the integration of temperature from composition
+    // is significantly more stable for exothermic systems
+    T += deltaT*R_[nSpecie];
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H
index 18e9745ac2..49541dd670 100644
--- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H
+++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,11 @@ Class
 Description
     An Euler implicit solver for chemistry
 
+    Euler implicit integration based on the reaction rate Jacobian is used to
+    solver for the composition and Euler explicit integration for the
+    temperature.  Separating the integration of temperature from composition
+    in this manner is significantly more stable for exothermic systems
+
 SourceFiles
     EulerImplicit.C
 
@@ -36,16 +41,13 @@ SourceFiles
 #define EulerImplicit_H
 
 #include "chemistrySolver.H"
-#include "Switch.H"
+#include "simpleMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-template <class Type>
-class simpleMatrix;
-
 /*---------------------------------------------------------------------------*\
                        Class EulerImplicit Declaration
 \*---------------------------------------------------------------------------*/
@@ -60,18 +62,21 @@ class EulerImplicit
         //- Coefficients dictionary
         dictionary coeffsDict_;
 
+        //- Chemistry timescale coefficient
+        scalar cTauChem_;
 
-        // Model constants
-
-            //- Chemistry timescale
-            scalar cTauChem_;
-
-            //- Equilibrium rate limiter flag (on/off)
-            Switch eqRateLimiter_;
-
-        // Solver data
+        //- Field encapsulating the composition, temperature and pressure
         mutable scalarField cTp_;
 
+        //- Reaction rate field
+        mutable scalarField R_;
+
+        //- Reaction Jacobian
+        mutable scalarSquareMatrix J_;
+
+        //- Euler implicit integration matrix for composition
+        mutable simpleMatrix<scalar> E_;
+
 
 public:
 
@@ -91,19 +96,6 @@ public:
 
     // Member Functions
 
-    void updateRRInReactionI
-    (
-        const label index,
-        const scalar pr,
-        const scalar pf,
-        const scalar corr,
-        const label lRef,
-        const label rRef,
-        const scalar p,
-        const scalar T,
-        simpleMatrix<scalar>& RR
-    ) const;
-
         //- Update the concentrations and return the chemical time
         virtual void solve
         (
diff --git a/tutorials/combustion/chemFoam/gri/constant/chemistryProperties b/tutorials/combustion/chemFoam/gri/constant/chemistryProperties
index a04784fbd1..742212d877 100644
--- a/tutorials/combustion/chemFoam/gri/constant/chemistryProperties
+++ b/tutorials/combustion/chemFoam/gri/constant/chemistryProperties
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-7;
 EulerImplicitCoeffs
 {
     cTauChem        1;
-    equilibriumRateLimiter off;
 }
 
 odeCoeffs
diff --git a/tutorials/combustion/chemFoam/h2/constant/chemistryProperties b/tutorials/combustion/chemFoam/h2/constant/chemistryProperties
index a04784fbd1..742212d877 100644
--- a/tutorials/combustion/chemFoam/h2/constant/chemistryProperties
+++ b/tutorials/combustion/chemFoam/h2/constant/chemistryProperties
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-7;
 EulerImplicitCoeffs
 {
     cTauChem        1;
-    equilibriumRateLimiter off;
 }
 
 odeCoeffs
diff --git a/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties b/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties
index a04784fbd1..742212d877 100644
--- a/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties
+++ b/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-7;
 EulerImplicitCoeffs
 {
     cTauChem        1;
-    equilibriumRateLimiter off;
 }
 
 odeCoeffs
diff --git a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
index 25b51b4583..52cc3c2a65 100644
--- a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
+++ b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-7;
 EulerImplicitCoeffs
 {
     cTauChem        10;
-    equilibriumRateLimiter no;
 }
 
 odeCoeffs
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
index 2e75e3f464..3b2ace0a6a 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
 EulerImplicitCoeffs
 {
     cTauChem        1;
-    equilibriumRateLimiter off;
 }
 
 odeCoeffs
diff --git a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
index 2e75e3f464..3b2ace0a6a 100644
--- a/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
+++ b/tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
 EulerImplicitCoeffs
 {
     cTauChem        1;
-    equilibriumRateLimiter off;
 }
 
 odeCoeffs
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/constant/gas/chemistryProperties b/tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/constant/gas/chemistryProperties
index a4201405a5..bece2af7b7 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/constant/gas/chemistryProperties
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/constant/gas/chemistryProperties
@@ -33,7 +33,6 @@ odeCoeffs
 EulerImplicitCoeffs
 {
     cTauChem        1;
-    equilibriumRateLimiter off;
 }
 
 #include "reactions"
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties
index c73670ef0d..6bae57ee20 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
 EulerImplicitCoeffs
 {
     cTauChem        0.05;
-    equilibriumRateLimiter off;
 }
 
 odeCoeffs
diff --git a/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
index 2e75e3f464..3b2ace0a6a 100644
--- a/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
+++ b/tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
 EulerImplicitCoeffs
 {
     cTauChem        1;
-    equilibriumRateLimiter off;
 }
 
 odeCoeffs
diff --git a/tutorials/multiphase/multiphaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas b/tutorials/multiphase/multiphaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
index b504cf02d3..709c77a8d5 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
+++ b/tutorials/multiphase/multiphaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
 EulerImplicitCoeffs
 {
     cTauChem        1;
-    equilibriumRateLimiter off;
 }
 
 odeCoeffs
diff --git a/tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/constant/chemistryProperties.vapor b/tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/constant/chemistryProperties.vapor
index d726e3d637..23647dd000 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/constant/chemistryProperties.vapor
+++ b/tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/constant/chemistryProperties.vapor
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
 EulerImplicitCoeffs
 {
     cTauChem                1;
-    equilibriumRateLimiter  off;
 }
 
 odeCoeffs
diff --git a/tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/constant/chemistryProperties.vapor b/tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/constant/chemistryProperties.vapor
index d726e3d637..23647dd000 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/constant/chemistryProperties.vapor
+++ b/tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/constant/chemistryProperties.vapor
@@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
 EulerImplicitCoeffs
 {
     cTauChem                1;
-    equilibriumRateLimiter  off;
 }
 
 odeCoeffs