diff --git a/include/cantera/numerics/DenseMatrix.h b/include/cantera/numerics/DenseMatrix.h
index 786cbae98..9861c8210 100644
--- a/include/cantera/numerics/DenseMatrix.h
+++ b/include/cantera/numerics/DenseMatrix.h
@@ -176,7 +176,7 @@ public:
 
     //  Listing of friend functions which are defined below
 
-    friend int solve(DenseMatrix& A, double* b);
+    friend int solve(DenseMatrix& A, double* b, size_t nrhs, size_t ldb);
     friend int solve(DenseMatrix& A, DenseMatrix& b);
     friend int invert(DenseMatrix& A, int nn);
 };
@@ -194,10 +194,12 @@ public:
  *
  *  The system is then solved using the LAPACK routine dgetrs
  *
- *   @param A   Dense matrix to be factored
- *   @param b   rhs to be solved.
+ *   @param A    Dense matrix to be factored
+ *   @param b    rhs(s) to be solved.
+ *   @param nrhs Number of right hand sides to solve
+ *   @param ldb  Leading dimension of b, if nrhs > 1
  */
-int solve(DenseMatrix& A, double* b);
+int solve(DenseMatrix& A, double* b, size_t nrhs=1, size_t ldb=0);
 
 //!  Solve Ax = b for multiple right-hand-side vectors.
 /*!
diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp
index 95462fe91..dd13191ed 100644
--- a/src/numerics/DenseMatrix.cpp
+++ b/src/numerics/DenseMatrix.cpp
@@ -131,7 +131,7 @@ vector_int& DenseMatrix::ipiv()
     return m_ipiv;
 }
 
-int solve(DenseMatrix& A, double* b)
+int solve(DenseMatrix& A, double* b, size_t nrhs, size_t ldb)
 {
     int info = 0;
     if (A.nColumns() != A.nRows()) {
@@ -164,8 +164,11 @@ int solve(DenseMatrix& A, double* b)
                             "DGETRF returned INFO = "+int2str(info) + ". The argument i has an illegal value");
     }
 
-    ct_dgetrs(ctlapack::NoTranspose, A.nRows(), 1, A.ptrColumn(0),
-              A.nRows(), &A.ipiv()[0], b, A.nColumns(), info);
+    if (ldb == 0) {
+        ldb = A.nColumns();
+    }
+    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) {
             writelogf("solve(DenseMatrix& A, double* b): DGETRS returned INFO = %d\n", info);
@@ -179,55 +182,7 @@ int solve(DenseMatrix& A, double* b)
 
 int solve(DenseMatrix& A, DenseMatrix& b)
 {
-    int info = 0;
-    if (A.nColumns() != A.nRows()) {
-        if (A.m_printLevel) {
-            writelogf("solve(DenseMatrix& A, DenseMatrix& b): Can only solve a square matrix\n");
-        }
-        if (! A.m_useReturnErrorCode) {
-            throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)", "Can only solve a square matrix");
-        }
-        return -1;
-    }
-    ct_dgetrf(A.nRows(), A.nColumns(), A.ptrColumn(0),
-              A.nRows(), &A.ipiv()[0], info);
-    if (info != 0) {
-        if (info > 0) {
-            if (A.m_printLevel) {
-                writelogf("solve(DenseMatrix& A, DenseMatrix& 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 CELapackError("solve(DenseMatrix& A, DenseMatrix& b)",
-                                    "DGETRF returned INFO = "+int2str(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.");
-            }
-        } else {
-            if (A.m_printLevel) {
-                writelogf("solve(DenseMatrix& A, DenseMatrix& b): DGETRF returned INFO = %d. The argument i has an illegal value\n", info);
-            }
-            if (! A.m_useReturnErrorCode) {
-                throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)",
-                                    "DGETRF returned INFO = "+int2str(info) + ". The argument i has an illegal value");
-            }
-        }
-        return info;
-    }
-
-    ct_dgetrs(ctlapack::NoTranspose, A.nRows(), b.nColumns(), A.ptrColumn(0),
-              A.nRows(), &A.ipiv()[0], b.ptrColumn(0), b.nRows(), info);
-    if (info != 0) {
-        if (A.m_printLevel) {
-            writelogf("solve(DenseMatrix& A, DenseMatrix& b): DGETRS returned INFO = %d\n", info);
-        }
-        if (! A.m_useReturnErrorCode) {
-            throw CELapackError("solve(DenseMatrix& A, DenseMatrix& b)", "DGETRS returned INFO = "+int2str(info));
-        }
-    }
-
-    return info;
+    return solve(A, b.ptrColumn(0), b.nColumns(), b.nRows());
 }
 
 void multiply(const DenseMatrix& A, const double* const b, double* const prod)