diff --git a/include/cantera/base/utilities.h b/include/cantera/base/utilities.h
index ff042173d..563444d89 100644
--- a/include/cantera/base/utilities.h
+++ b/include/cantera/base/utilities.h
@@ -521,6 +521,15 @@ void deepStdVectorPointerCopy(const std::vector<D*> &fromVec, std::vector<D*> &t
 //! Check to see that a number is finite (not NaN, +Inf or -Inf)
 void checkFinite(const double tmp);
 
+//! Check to see that all elements in an array are finite
+/*!
+ *  Throws an exception if any element is NaN, +Inf, or -Inf
+ *  @param name    Name to be used in the exception message if the check fails
+ *  @param values  Array of *N* values to be checked
+ *  @param N       Number of elements in *values*
+ */
+void checkFinite(const std::string& name, double* values, size_t N);
+
 //! Const accessor for a value in a std::map.
 /*!
  *  This is a const alternative to operator[]. Roughly equivalent to the 'at'
diff --git a/src/base/checkFinite.cpp b/src/base/checkFinite.cpp
index 316e22757..134679eaf 100644
--- a/src/base/checkFinite.cpp
+++ b/src/base/checkFinite.cpp
@@ -13,6 +13,8 @@
 
 #include <stdexcept>
 #include <cstdio>
+#include "cantera/base/stringUtils.h"
+#include "cantera/base/ctexceptions.h"
 
 // We expect that there will be special casing based on the computer
 // system here
@@ -54,4 +56,20 @@ void checkFinite(const double tmp)
     }
 }
 
+void checkFinite(const std::string& name, double* values, size_t N)
+{
+    for (size_t i = 0; i < N; i++) {
+        if (!finite(values[i])) {
+            std::string message = name + " contains non-finite elements:\n\n";
+            for (size_t j = 0; j < N; j++) {
+                if (!finite(values[j])) {
+                    message += name + "[" + int2str(j) + "] = " +
+                               fp2str(values[j]) + "\n";
+                }
+            }
+            throw CanteraError("checkFinite", message);
+        }
+    }
+}
+
 }
diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp
index 560a244c4..88852372b 100644
--- a/src/zeroD/IdealGasReactor.cpp
+++ b/src/zeroD/IdealGasReactor.cpp
@@ -56,11 +56,6 @@ void IdealGasReactor::initialize(doublereal t0)
 
 void IdealGasReactor::updateState(doublereal* y)
 {
-    for (size_t i = 0; i < m_nv; i++) {
-        AssertFinite(y[i], "IdealGasReactor::updateState",
-                     "y[" + int2str(i) + "] is not finite");
-    }
-
     // The components of y are [0] the total mass, [1] the total volume,
     // [2] the temperature, [3...K+3] are the mass fractions of each species,
     // and [K+3...] are the coverages of surface species on each wall.
@@ -142,11 +137,6 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y,
         ydot[2] = 0;
     }
 
-    for (size_t i = 0; i < m_nv; i++) {
-        AssertFinite(ydot[i], "IdealGasReactor::evalEqs",
-                     "ydot[" + int2str(i) + "] is not finite");
-    }
-
     resetSensitivity(params);
 }
 
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index 5e3f0a145..b47433ff9 100644
--- a/src/zeroD/Reactor.cpp
+++ b/src/zeroD/Reactor.cpp
@@ -128,11 +128,6 @@ void Reactor::syncState()
 
 void Reactor::updateState(doublereal* y)
 {
-    for (size_t i = 0; i < m_nv; i++) {
-        AssertFinite(y[i], "Reactor::updateState",
-                     "y[" + int2str(i) + "] is not finite");
-    }
-
     // The components of y are [0] the total mass, [1] the total volume,
     // [2] the total internal energy, [3...K+3] are the mass fractions of each
     // species, and [K+3...] are the coverages of surface species on each wall.
@@ -267,11 +262,6 @@ void Reactor::evalEqs(doublereal time, doublereal* y,
     }
 
     ydot[0] = dmdt;
-
-    for (size_t i = 0; i < m_nv; i++) {
-        AssertFinite(ydot[i], "Reactor::evalEqs",
-                     "ydot[" + int2str(i) + "] is not finite");
-    }
     resetSensitivity(params);
 }
 
diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp
index 8f429dae6..490fe9ce8 100644
--- a/src/zeroD/ReactorNet.cpp
+++ b/src/zeroD/ReactorNet.cpp
@@ -156,6 +156,7 @@ void ReactorNet::eval(doublereal t, doublereal* y,
                                ydot + m_start[n], p + pstart);
         pstart += m_nparams[n];
     }
+    checkFinite("ydot", ydot, m_nv);
 }
 
 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
@@ -186,6 +187,7 @@ void ReactorNet::evalJacobian(doublereal t, doublereal* y,
 
 void ReactorNet::updateState(doublereal* y)
 {
+    checkFinite("y", y, m_nv);
     for (size_t n = 0; n < m_reactors.size(); n++) {
         m_reactors[n]->updateState(y + m_start[n]);
     }