diff --git a/include/cantera/numerics/BandMatrix.h b/include/cantera/numerics/BandMatrix.h
index a7e8a4f5c..1c04dfd03 100644
--- a/include/cantera/numerics/BandMatrix.h
+++ b/include/cantera/numerics/BandMatrix.h
@@ -280,6 +280,9 @@ public:
      */
     virtual size_t checkColumns(doublereal& valueSmall) const;
 
+    //! LAPACK "info" flag after last factor/solve operation
+    int info() const { return m_info; };
+
 protected:
     //! Matrix data
     vector_fp data;
@@ -311,6 +314,8 @@ protected:
 
     //! Extra dp work array needed - size = 3n
     vector_fp work_;
+
+    int m_info;
 };
 
 //! Utility routine to print out the matrix
diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp
index 4830d257b..c2ec04e7d 100644
--- a/src/numerics/BandMatrix.cpp
+++ b/src/numerics/BandMatrix.cpp
@@ -30,7 +30,8 @@ BandMatrix::BandMatrix() :
     m_n(0),
     m_kl(0),
     m_ku(0),
-    m_zero(0.0)
+    m_zero(0.0),
+    m_info(0)
 {
 }
 
@@ -38,7 +39,8 @@ BandMatrix::BandMatrix(size_t n, size_t kl, size_t ku, doublereal v)   :
     m_n(n),
     m_kl(kl),
     m_ku(ku),
-    m_zero(0.0)
+    m_zero(0.0),
+    m_info(0)
 {
     data.resize(n*(2*kl + ku + 1));
     ludata.resize(n*(2*kl + ku + 1));
@@ -59,7 +61,8 @@ BandMatrix::BandMatrix(const BandMatrix& y) :
     m_n(0),
     m_kl(0),
     m_ku(0),
-    m_zero(0.0)
+    m_zero(0.0),
+    m_info(y.m_info)
 {
     m_n = y.m_n;
     m_kl = y.m_kl;
@@ -95,6 +98,7 @@ BandMatrix& BandMatrix::operator=(const BandMatrix& y)
         m_colPtrs[j] = &data[ldab * j];
         m_lu_col_ptrs[j] = &ludata[ldab * j];
     }
+    m_info = y.m_info;
     return *this;
 }
 
@@ -235,27 +239,23 @@ void BandMatrix::leftMult(const doublereal* const b, doublereal* const prod) con
 
 int BandMatrix::factor()
 {
-    int info=0;
     ludata = data;
 #if CT_USE_LAPACK
     ct_dgbtrf(nRows(), nColumns(), nSubDiagonals(), nSuperDiagonals(),
-              ludata.data(), ldim(), ipiv().data(), info);
+              ludata.data(), ldim(), ipiv().data(), m_info);
 #else
     long int nu = static_cast<long int>(nSuperDiagonals());
     long int nl = static_cast<long int>(nSubDiagonals());
     long int smu = nu + nl;
-    info = bandGBTRF(m_lu_col_ptrs.data(), static_cast<long int>(nColumns()),
-                     nu, nl, smu, m_ipiv.data());
+    m_info = bandGBTRF(m_lu_col_ptrs.data(), static_cast<long int>(nColumns()),
+                       nu, nl, smu, m_ipiv.data());
 #endif
-    // if info = 0, LU decomp succeeded.
-    if (info == 0) {
-        m_factored = true;
-    } else {
-        m_factored = false;
-        ofstream fout("bandmatrix.csv");
-        fout << *this << endl;
+    if (m_info != 0) {
+        throw Cantera::CanteraError("BandMatrix::factor",
+            "Factorization failed with DGBTRF error code {}.", m_info);
     }
-    return info;
+    m_factored = true;
+    return m_info;
 }
 
 int BandMatrix::solve(const doublereal* const b, doublereal* const x)
@@ -266,34 +266,30 @@ int BandMatrix::solve(const doublereal* const b, doublereal* const x)
 
 int BandMatrix::solve(doublereal* b, size_t nrhs, size_t ldb)
 {
-    int info = 0;
     if (!m_factored) {
-        info = factor();
+        factor();
     }
     if (ldb == 0) {
         ldb = nColumns();
     }
-    if (info == 0) {
 #if CT_USE_LAPACK
-        ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(),
-                  nSuperDiagonals(), nrhs, ludata.data(), ldim(),
-                  ipiv().data(), b, ldb, info);
+    ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(),
+              nSuperDiagonals(), nrhs, ludata.data(), ldim(),
+              ipiv().data(), b, ldb, m_info);
 #else
-        long int nu = static_cast<long int>(nSuperDiagonals());
-        long int nl = static_cast<long int>(nSubDiagonals());
-        long int smu = nu + nl;
-        double** a = m_lu_col_ptrs.data();
-        bandGBTRS(a, static_cast<long int>(nColumns()), smu, nl, m_ipiv.data(),
-                  b);
+    long int nu = static_cast<long int>(nSuperDiagonals());
+    long int nl = static_cast<long int>(nSubDiagonals());
+    long int smu = nu + nl;
+    double** a = m_lu_col_ptrs.data();
+    bandGBTRS(a, static_cast<long int>(nColumns()), smu, nl, m_ipiv.data(), b);
+    m_info = 0;
 #endif
-    }
 
-    // error handling
-    if (info != 0) {
-        ofstream fout("bandmatrix.csv");
-        fout << *this << endl;
+    if (m_info != 0) {
+        throw Cantera::CanteraError("BandMatrix::solve",
+            "Linear solve failed with DGBTRS error code {}.", m_info);
     }
-    return info;
+    return m_info;
 }
 
 vector_fp::iterator BandMatrix::begin()
diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp
index fbb07c0ab..ca97e8820 100644
--- a/src/oneD/MultiNewton.cpp
+++ b/src/oneD/MultiNewton.cpp
@@ -167,27 +167,29 @@ void MultiNewton::step(doublereal* x, doublereal* step,
         step[n] = -step[n];
     }
 
-    size_t iok = jac.solve(step, step);
-    // if iok is non-zero, then solve failed
-    if (iok != 0) {
-        iok--;
-        size_t nd = r.nDomains();
-        size_t n;
-        for (n = nd-1; n != npos; n--) {
-            if (iok >= r.start(n)) {
-                break;
+    try {
+        jac.solve(step, step);
+    } catch (CanteraError&) {
+        int iok = jac.info() - 1;
+        if (iok >= 0) {
+            size_t nd = r.nDomains();
+            size_t n;
+            for (n = nd-1; n != npos; n--) {
+                if (iok >= static_cast<int>(r.start(n))) {
+                    break;
+                }
             }
+            Domain1D& dom = r.domain(n);
+            size_t offset = iok - r.start(n);
+            size_t pt = offset/dom.nComponents();
+            size_t comp = offset - pt*dom.nComponents();
+            throw CanteraError("MultiNewton::step",
+                "Jacobian is singular for domain {}, component {} at point {}\n"
+                "(Matrix row {})",
+                dom.id(), dom.componentName(comp), pt, iok);
+        } else {
+            throw;
         }
-        Domain1D& dom = r.domain(n);
-        size_t offset = iok - r.start(n);
-        size_t pt = offset/dom.nComponents();
-        size_t comp = offset - pt*dom.nComponents();
-        throw CanteraError("MultiNewton::step",
-            "Jacobian is singular for domain {}, component {} at point {}\n"
-            "(Matrix row {}) \nsee file bandmatrix.csv\n",
-            dom.id(), dom.componentName(comp), pt, iok);
-    } else if (int(iok) < 0) {
-        throw CanteraError("MultiNewton::step", "iok = {}", iok);
     }
 }