diff --git a/include/cantera/numerics/DenseMatrix.h b/include/cantera/numerics/DenseMatrix.h
index 409041ce33..32afda839a 100644
--- a/include/cantera/numerics/DenseMatrix.h
+++ b/include/cantera/numerics/DenseMatrix.h
@@ -155,6 +155,8 @@ class DenseMatrix : public Array2D
     friend int invert(DenseMatrix& A, int nn);
 };
 
+int factor(DenseMatrix& A);
+int solveFactored(DenseMatrix& A, double* b, size_t nrhs=1, size_t ldb=0);
 
 //! Solve Ax = b. Array b is overwritten on exit with x.
 /*!
diff --git a/include/cantera/numerics/Newton.h b/include/cantera/numerics/Newton.h
index d68c23f60c..b827937070 100644
--- a/include/cantera/numerics/Newton.h
+++ b/include/cantera/numerics/Newton.h
@@ -52,6 +52,10 @@ class Newton
     //! Compute the weighted 2-norm of `step`.
     doublereal weightedNorm(const doublereal* x, const doublereal* step) const;
 
+    int hybridSolve();
+
+    int timestep();
+
     /**
      * Find the solution to F(X) = 0 by damped Newton iteration.
      */
@@ -75,32 +79,46 @@ class Newton
         m_max[n] = upper;
     }
 
-    void setConstant(size_t component, bool constant) {
-        m_constant[component] = constant;
+    void setConstants(vector_int constantComponents) {
+        m_constantComponents = constantComponents;
     }
 
     void evalJacobian(doublereal* x, doublereal* xdot);
 
+    void getSolution(double* x) {
+        for (size_t i = 0; i < m_nv; i++) {
+            x[i] = m_x[i];
+        }
+    }
+
 protected:
     FuncEval* m_residfunc;
 
-    DenseMatrix m_jacobian;
+    //! number of variables
+    size_t m_nv;
+
+    //! solution converged if [weightedNorm(sol, step) < m_convergenceThreshold]
+    doublereal m_convergenceThreshold;
+
+    DenseMatrix m_jacobian, m_jacFactored;
     int m_jacAge, m_jacMaxAge;
     doublereal m_jacRtol, m_jacAtol;
 
 
-    //! Work arrays of size #m_nv used in solve().
+    //! work arrays of size #m_nv used in solve().
     vector_fp m_x, m_x1, m_stp, m_stp1;
 
     vector_fp m_max, m_min;
     vector_fp m_rtol_ss, m_rtol_ts;
     vector_fp m_atol_ss, m_atol_ts;
 
-    //! number of variables
-    size_t m_nv;
+    vector_fp m_xlast, m_xsave;
+
+    //! the indexes of any constant variables
+    vector_int m_constantComponents;
 
-    //! constant variables
-    std::vector<bool> m_constant;
+    //! current timestep reciprocal
+    doublereal m_rdt = 0;
 };
 }
 
diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp
index 2d5b014e01..e9d6783cec 100644
--- a/src/numerics/DenseMatrix.cpp
+++ b/src/numerics/DenseMatrix.cpp
@@ -140,6 +140,83 @@ vector_int& DenseMatrix::ipiv()
     return m_ipiv;
 }
 
+int factor(DenseMatrix& A)
+{
+    int info = 0;
+    #if CT_USE_LAPACK
+        ct_dgetrf(A.nRows(), A.nColumns(), A.ptrColumn(0),
+                  A.nRows(), &A.ipiv()[0], info);
+        if (info > 0) {
+            if (A.m_printLevel) {
+                writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d   U(i,i) is exactly zero. The factorization has"
+                          " been completed, but the factor U is exactly singular, and division by zero will occur if "
+                          "it is used to solve a system of equations.\n", info);
+            }
+            if (!A.m_useReturnErrorCode) {
+                throw CanteraError("solve(DenseMatrix& A, double* b)",
+                                    "DGETRF returned INFO = {}. U(i,i) is exactly zero. The factorization has"
+                                    " been completed, but the factor U is exactly singular, and division by zero will occur if "
+                                    "it is used to solve a system of equations.", info);
+            }
+            return info;
+        } else if (info < 0) {
+            if (A.m_printLevel) {
+                writelogf("solve(DenseMatrix& A, double* b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info);
+            }
+
+            throw CanteraError("solve(DenseMatrix& A, double* b)",
+                               "DGETRF returned INFO = {}. The argument i has an illegal value", info);
+        }
+    #else
+        MappedMatrix Am(&A(0,0), A.nRows(), A.nColumns());
+        #ifdef NDEBUG
+            auto lu = Am.partialPivLu();
+        #else
+            auto lu = Am.fullPivLu();
+            if (lu.nonzeroPivots() < static_cast<long int>(A.nColumns())) {
+                throw CanteraError("solve(DenseMatrix& A, double* b)",
+                    "Matrix appears to be rank-deficient: non-zero pivots = {}; columns = {}",
+                    lu.nonzeroPivots(), A.nColumns());
+            }
+        #endif
+    #endif
+    return info;
+}
+
+int solveFactored(DenseMatrix& A, double* b, size_t nrhs, size_t ldb)
+{
+    if (A.nColumns() != A.nRows()) {
+        if (A.m_printLevel) {
+            writelogf("solve(DenseMatrix& A, double* b): Can only solve a square matrix\n");
+        }
+        throw CanteraError("solve(DenseMatrix& A, double* b)", "Can only solve a square matrix");
+    }
+
+    int info = 0;
+    if (ldb == 0) {
+        ldb = A.nColumns();
+    }
+    #if CT_USE_LAPACK
+        ct_dgetrs(ctlapack::NoTranspose, A.nRows(), nrhs, A.ptrColumn(0),
+                  A.nRows(), &A.ipiv()[0], b, ldb, info);
+        if (info != 0) {
+            if (A.m_printLevel) {
+                writelog("solve(DenseMatrix& A, double* b): DGETRS returned INFO = {}\n", info);
+            }
+            if (info < 0 || !A.m_useReturnErrorCode) {
+                throw CanteraError("solve(DenseMatrix& A, double* b)",
+                                   "DGETRS returned INFO = {}", info);
+            }
+        }
+    #else
+        for (size_t i = 0; i < nrhs; i++) {
+            MappedVector bm(b + ldb*i, A.nColumns());
+            bm = lu.solve(bm);
+        }
+    #endif
+    return info;
+}
+
 int solve(DenseMatrix& A, double* b, size_t nrhs, size_t ldb)
 {
     if (A.nColumns() != A.nRows()) {
diff --git a/src/numerics/Newton.cpp b/src/numerics/Newton.cpp
index b5eab8f7ec..cdbb5917f6 100644
--- a/src/numerics/Newton.cpp
+++ b/src/numerics/Newton.cpp
@@ -20,7 +20,7 @@ const size_t NDAMP = 7;
 Newton::Newton(FuncEval& func) {
     m_residfunc = &func;
     m_nv = m_residfunc->neq();
-    m_constant.resize(m_nv, false);
+    m_convergenceThreshold = 1.0e-14;
     m_x.resize(m_nv);
     m_x1.resize(m_nv);
     m_stp.resize(m_nv);
@@ -31,93 +31,102 @@ Newton::Newton(FuncEval& func) {
     m_atol_ss.resize(m_nv, 1.0e-9);
     // m_rtol_ts.resize(m_nv, 1.0e-4);
     // m_atol_ts.resize(m_nv, 1.0e-11);
+    m_xlast.resize(m_nv);
+    m_xsave.resize(m_nv);
+
+    m_rdt = 0;
 
     m_jacobian = DenseMatrix(m_nv, m_nv);
     m_jacAge = 10000;
-    m_jacMaxAge = 5;
-    m_jacRtol = 1.0e-5;
+    m_jacMaxAge = 1;
+    m_jacRtol = 1.0e-15;
     m_jacAtol = sqrt(std::numeric_limits<double>::epsilon());
 }
 
 void Newton::evalJacobian(doublereal* x, doublereal* xdot) {
 
-    // // calculate unperturbed residual
+    // calculate unperturbed residual
     m_residfunc->eval(0, x, xdot, 0);
 
     for (size_t n = 0; n < m_nv; n++) {
-        // calculate the nth Jacobian column, unless component n is constant
-        if (!m_constant[n]) {
-            double xsave = x[n];
-
-            // calculate the perturbation amount, preserving the sign of x[n]
-            double dx;
-            if (xsave >= 0) {
-                dx = xsave*m_jacRtol + m_jacAtol;
-            } else {
-                dx = xsave*m_jacRtol - m_jacAtol;
-            }
+        // calculate the nth Jacobian column
+        double xsave = x[n];
 
-            // perturb the solution vector
-            x[n] = xsave + dx;
-            dx = x[n] - xsave;
+        // calculate the perturbation amount, preserving the sign of x[n]
+        double dx;
+        if (xsave >= 0) {
+            dx = xsave*m_jacRtol + m_jacAtol;
+        } else {
+            dx = xsave*m_jacRtol - m_jacAtol;
+        }
 
-            // calculate perturbed residual
-            vector_fp xdotPerturbed(m_nv); //make this member for speed?
-            m_residfunc->eval(0, x, xdotPerturbed.data(), 0);
+        // perturb the solution vector
+        x[n] = xsave + dx;
+        dx = x[n] - xsave;
 
-            // compute nth column of Jacobian
-            for (size_t m = 0; m < m_nv; m++) {
-                m_jacobian.value(m,n) = (xdotPerturbed[m] - xdot[m])/dx;
-            }
-            // restore solution vector
-            x[n] = xsave;
+        // calculate perturbed residual
+        vector_fp xdotPerturbed(m_nv); //make this member for speed?
+        m_residfunc->eval(0, x, xdotPerturbed.data(), 0);
 
-        // for constant components: Jacobian column of 0's, with 1 on diagonal
-        // note: Jacobian row must also be 0's w/ 1 on diagonal
-        } else {
-            for (size_t m = 0; m < m_nv; m++) {
-                m_jacobian.value(m,n) = 0;
+        // compute nth column of Jacobian
+        for (size_t m = 0; m < m_nv; m++) {
+            m_jacobian.value(m,n) = (xdotPerturbed[m] - xdot[m])/dx;
+        }
+        // restore solution vector
+        x[n] = xsave;
+
+        // for constant-valued components: 1 in diagonal position, all 0's in row and column
+        for (size_t i : m_constantComponents) {
+            for (size_t j = 0; j < m_nv; j++) {
+                m_jacobian.value(i,j) = 0;
+                m_jacobian.value(j,i) = 0;
             }
-            m_jacobian.value(n,n) = 1;
+            m_jacobian.value(i,i) = 1;
         }
     }
 
     // writelog("\nnew jac:\n");
     // for (int i = 0; i < m_nv; i++) {
+    //     writelog("ROW {} | ", i);
     //     for (int j = 0; j < m_nv; j++) {
-    //         writelog("{:14.5} ", m_jacobian.value(i,j));
+    //         writelog("{:.5} ", m_jacobian.value(i,j));
     //     }
-    //     writelog("\n");
+    //     writelog("\n\n");
     // }
 
     m_jacAge = 0;
+
+    m_jacFactored = m_jacobian;
+    Cantera::factor(m_jacFactored);
 }
 
 doublereal Newton::weightedNorm(const doublereal* x, const doublereal* step) const
 {
-    double sum = 0.0;
+    double accum = 0.0;
     for (size_t n = 0; n < m_nv; n++) {
-        double weight = m_rtol_ss[n]*fabs(x[n]) + m_atol_ss[n]; //TODO: add transient tolerances if rdt!=0
-        double f = step[n]/weight;
-        sum += f*f;
+        double f = step[n]/(x[n] + m_atol_ss[n]);
+        accum += f*f;
     }
-    return sqrt(sum/m_nv);
+    // writelog("test.... {}", m_constantComponents.size());
+    // return sqrt(accum/(m_nv - m_constantComponents.size()));
+    return sqrt(accum/m_nv);
 }
 
 void Newton::step(doublereal* x, doublereal* step, int loglevel)
 {
     m_residfunc->eval(0, x, step, 0);
     for (size_t n = 0; n < m_nv; n++) {
-        step[n] = -step[n];
+        step[n] = -step[n] + m_rdt*(x[n]-m_xlast[n]);
     }
 
-    DenseMatrix solvejac = m_jacobian;
+    DenseMatrix solvejac = m_jacFactored;
     try {
         //Note: this function takes an unfactored jacobian, then finds its LU factorization before
         // solving. Optimization is possible by saving the factored jacobian, since it can be reused.
         // Also, the DenseMatrix provided here will be overwritten with the LU factored version, so
         // a copy is passed instead in order to preserve the original for reuse.
-        Cantera::solve(solvejac, step);
+        // Cantera::solve(solvejac, step);
+        Cantera::solveFactored(solvejac, step);
 
     } catch (CanteraError&) {
         // int iok = m_jac->info() - 1;
@@ -222,9 +231,9 @@ int Newton::dampStep(const doublereal* x0, const doublereal* step0,
         // write log information
         if (loglevel > 0) {
             doublereal ss = weightedNorm(x1,step1);
-            writelog("\n{:d}  {:9.5f}   {:9.5f}   {:9.5f}   {:9.5f}   {:9.5f}  {:d}/{:d}",
-                     m, damp, boundFactor, log10(ss+SmallNumber),
-                     log10(s0+SmallNumber), log10(s1+SmallNumber),
+            writelog("\n{:d}  {:9.5f}   {:9.5f}   {:9.5e}   {:9.5e}   {:9.5e}  {:d}/{:d}",
+                     m, damp, boundFactor, ss+SmallNumber,
+                     s0+SmallNumber, s1+SmallNumber,
                      m_jacAge, m_jacMaxAge);
         }
 
@@ -232,7 +241,7 @@ int Newton::dampStep(const doublereal* x0, const doublereal* step0,
         // damping coefficient. Also accept it if this step would result in a
         // converged solution. Otherwise, decrease the damping coefficient and
         // try again.
-        if (s1 < 1.0 || s1 < s0) {
+        if (s1 < m_convergenceThreshold || s1 < s0) {
             break;
         }
         damp /= DampFactor;
@@ -242,7 +251,7 @@ int Newton::dampStep(const doublereal* x0, const doublereal* step0,
     // stepping by the damped step would represent a converged solution, and
     // return 0 otherwise. If no damping coefficient could be found, return -2.
     if (m < NDAMP) {
-        if (s1 > 1.0) {
+        if (s1 > m_convergenceThreshold) {
             return 0;
         } else {
             return 1;
@@ -252,6 +261,44 @@ int Newton::dampStep(const doublereal* x0, const doublereal* step0,
     }
 }
 
+int Newton::hybridSolve() {
+    int MAX = 100;
+    int newtonsolves = 0;
+    int timesteps = 0;
+
+    m_residfunc->getState(m_x.data());
+    copy(m_x.begin(), m_x.end(), m_xsave.begin());
+
+    for(int i = 0; i < MAX; i++) {
+        newtonsolves++;
+        if(solve() == 1) {
+            writelog("\nConverged in {} newton solves, {} timesteps.", newtonsolves, timesteps);
+            return 1;
+        }
+        copy(m_xsave.begin(), m_xsave.end(), m_x.begin());
+        for(int j = 0; j < MAX; j++) {
+            timesteps++;
+            timestep();
+        }
+        copy(m_x.begin(), m_x.end(), m_xsave.begin());
+    }
+    writelog("Solver failure...");
+    return 0;
+}
+
+int Newton::timestep() {
+    m_rdt = 1.0/1.0e-5;
+
+    // calculate time-integration Jacobian
+    // m_residfunc->getState(m_x.data());
+    copy(m_x.begin(), m_x.end(), m_xlast.begin());
+    evalJacobian(&m_x[0], &m_stp[0]);
+    for (size_t i = 0; i < m_nv; i++) {
+        m_jacobian.value(i,i) -= m_rdt;
+    }
+    return solve();
+}
+
 int Newton::solve(int loglevel)
 {
     int m = 0;
@@ -302,6 +349,16 @@ int Newton::solve(int loglevel)
         if (m == 0) {
             copy(m_x1.begin(), m_x1.end(), m_x.begin());
         } else if (m == 1) {
+            // writelog("\nConverged. Newton steps: {}", steps);
+            // writelog("\nconverged!\n");
+            // writelog("\nsolution components: ");
+            // for(int i = 0; i < m_nv; i++) {
+            //     writelog("{:9.5e}, ", m_x[i]);
+            // }
+            // writelog("\nresidual components: ");
+            // for(int i = 0; i < m_nv; i++) {
+            //     writelog("{:9.5e}, ", m_stp[i]);
+            // }
             // convergence
             // if (rdt == 0) {
             //     jac.setAge(0); // for efficient sensitivity analysis
diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp
index f032e7626d..28b44c91f8 100644
--- a/src/zeroD/ReactorNet.cpp
+++ b/src/zeroD/ReactorNet.cpp
@@ -419,13 +419,19 @@ double ReactorNet::solveSteady()
     {
         m_newt->setBounds(i, -1.0e-3, 1.01);
     }
-    m_newt->setBounds(0, 0, 100);
+    m_newt->setBounds(0, 1.0e-3, 100);
     m_newt->setBounds(1, 0, 5);
     m_newt->setBounds(2, -10000000, 10000000);
 
-    m_newt->setConstant(1, true);
+    m_newt->setConstants({1});
+    // m_newt->setConstants({0,1,2,11,12});
+    // m_newt->setConstant(0, true);
+    // m_newt->setConstant(1, true);
+    // m_newt->setConstant(2, true);
+    // m_newt->setConstant(11, true);
+    // m_newt->setConstant(12, true);
 
-    m_newt->solve(8);
+    return m_newt->hybridSolve();
 }
 
 }