diff --git a/include/cantera/numerics/BandMatrix.h b/include/cantera/numerics/BandMatrix.h
index c3a7af074..c7d21e454 100644
--- a/include/cantera/numerics/BandMatrix.h
+++ b/include/cantera/numerics/BandMatrix.h
@@ -191,7 +191,7 @@ public:
      *          0 indicates a success
      *         ~0  Some error occurred, see the LAPACK documentation
      */
-    int solve(doublereal* b);
+    int solve(doublereal* b, size_t nrhs=1, size_t ldb=0);
 
     //! Returns an iterator for the start of the band storage data
     /*!
diff --git a/include/cantera/numerics/GeneralMatrix.h b/include/cantera/numerics/GeneralMatrix.h
index 6e9374580..c59fcde46 100644
--- a/include/cantera/numerics/GeneralMatrix.h
+++ b/include/cantera/numerics/GeneralMatrix.h
@@ -127,9 +127,12 @@ public:
 
     //! Solves the Ax = b system returning x in the b spot.
     /*!
-     *  @param b  Vector for the rhs of the equation system
+     *  @param b    Vector for the rhs of the equation system
+     *  @param nrhs Number of right-hand sides to solve, default 1
+     *  @param ldb  Leading dimension of the right-hand side array.
+     *              Defaults to nRows()
      */
-    virtual int solve(doublereal* b) = 0;
+    virtual int solve(doublereal* b, size_t nrhs=1, size_t ldb=0) = 0;
 
     //! true if the current factorization is up to date with the matrix
     virtual bool factored() const = 0;
diff --git a/include/cantera/numerics/SquareMatrix.h b/include/cantera/numerics/SquareMatrix.h
index f71a37b1b..fc7133c04 100644
--- a/include/cantera/numerics/SquareMatrix.h
+++ b/include/cantera/numerics/SquareMatrix.h
@@ -44,7 +44,7 @@ public:
     //! Assignment operator
     SquareMatrix& operator=(const SquareMatrix& right);
 
-    int solve(doublereal* b);
+    int solve(doublereal* b, size_t nrhs=1, size_t ldb=0);
 
     void resize(size_t n, size_t m, doublereal v = 0.0);
 
diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp
index a449a9049..c7dc96be3 100644
--- a/src/numerics/BandMatrix.cpp
+++ b/src/numerics/BandMatrix.cpp
@@ -270,16 +270,19 @@ int BandMatrix::solve(const doublereal* const b, doublereal* const x)
     return solve(x);
 }
 
-int BandMatrix::solve(doublereal* b)
+int BandMatrix::solve(doublereal* b, size_t nrhs, size_t ldb)
 {
     int info = 0;
     if (!m_factored) {
         info = factor();
     }
+    if (ldb == 0) {
+        ldb = nColumns();
+    }
     if (info == 0)
         ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(),
-                  nSuperDiagonals(), 1, DATA_PTR(ludata), ldim(),
-                  DATA_PTR(ipiv()), b, nColumns(), info);
+                  nSuperDiagonals(), nrhs, DATA_PTR(ludata), ldim(),
+                  DATA_PTR(ipiv()), b, ldb, info);
 
     // error handling
     if (info != 0) {
diff --git a/src/numerics/SquareMatrix.cpp b/src/numerics/SquareMatrix.cpp
index 4e154e264..9959515b9 100644
--- a/src/numerics/SquareMatrix.cpp
+++ b/src/numerics/SquareMatrix.cpp
@@ -60,7 +60,7 @@ SquareMatrix& SquareMatrix::operator=(const SquareMatrix& y)
     return *this;
 }
 
-int SquareMatrix::solve(doublereal* b)
+int SquareMatrix::solve(doublereal* b, size_t nrhs, size_t ldb)
 {
     if (useQR_) {
         return solveQR(b);
@@ -75,12 +75,15 @@ int SquareMatrix::solve(doublereal* b)
             return retn;
         }
     }
+    if (ldb == 0) {
+        ldb = nColumns();
+    }
     /*
      * Solve the factored system
      */
     ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(nRows()),
-              1, &(*(begin())), static_cast<int>(nRows()),
-              DATA_PTR(ipiv()), b, static_cast<int>(nColumns()), info);
+              nrhs, &(*(begin())), static_cast<int>(nRows()),
+              DATA_PTR(ipiv()), b, ldb, info);
     if (info != 0) {
         if (m_printLevel) {
             writelogf("SquareMatrix::solve(): DGETRS returned INFO = %d\n", info);