diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 019fbaab3966006377af500645cd28ce1a163797..c2c158a2d8ef1693959aa59befd6aa56a7a35f32 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -847,7 +847,7 @@ auto trace_VLE_isotherm_binary(const Model &model, Scalar T, VecType rhovecL0, V
             auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[0]) + N, N);
             auto x = rhovecL / rhovecL.sum();
             auto y = rhovecV / rhovecV.sum();
-            if ((x < 0).any() || (x > 1).any() || (y < 0).any() || (y > 1).any()) {
+            if ((x < 0).any() || (x > 1).any() || (y < 0).any() || (y > 1).any() || (!rhovecL.isFinite()).any() || (!rhovecV.isFinite()).any()) {
                 return true;
             }
             else {
@@ -1076,7 +1076,7 @@ auto trace_VLE_isobar_binary(const Model& model, Scalar p, Scalar T0, VecType rh
             auto rhovecV = Eigen::Map<const Eigen::ArrayXd>(&(x0[1]) + N, N);
             auto x = rhovecL / rhovecL.sum();
             auto y = rhovecV / rhovecV.sum();
-            if ((x < 0).any() || (x > 1).any() || (y < 0).any() || (y > 1).any()) {
+            if ((x < 0).any() || (x > 1).any() || (y < 0).any() || (y > 1).any() || (!rhovecL.isFinite()).any() || (!rhovecV.isFinite()).any()) {
                 return true;
             }
             else {