diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index 0f97b368bd95df8aee92baec96a39f4f89f8ece2..99b9f355eddc012e0f85604b0fa35d7f93d44edf 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -443,6 +443,18 @@ struct CriticalTracing {
 
             double dtold = dt;
             auto x0_previous = x0;
+
+            // Call the function to get drho/dt in order 
+            // to cache it at the *beginning* of the step
+            auto dxdt = x0;
+            try {
+                xprime(x0, dxdt, -1.0);
+            }
+            catch (...) {
+                break;
+            }
+            auto drhodt = Eigen::Map<Eigen::ArrayXd>(&(dxdt[0]) + 1, dxdt.size() - 1);
+            last_drhodt = drhodt;
             
             if (options.integration_order == 5) {
                 controlled_step_result res = controlled_step_result::fail;
@@ -466,20 +478,17 @@ struct CriticalTracing {
                 dt = std::min(dt, options.max_dt);
             }
             else if (options.integration_order == 1) {
-                eul.do_step(xprime, x0, t, dt);
+                try {
+                    eul.do_step(xprime, x0, t, dt);
+                }
+                catch (...) {
+                    break;
+                }
             }
             else {
                 throw std::invalid_argument("integration order is invalid:" + std::to_string(options.integration_order));
             }
 
-            // Call the function again to get drho/dt in order 
-            // to cache it
-            auto dxdt = x0;
-            xprime(x0, dxdt, -1.0);
-            auto drhodt = Eigen::Map<Eigen::ArrayXd>(&(dxdt[0]) + 1, dxdt.size() - 1);
-            last_drhodt = drhodt;
-            auto absrelchange = std::abs(drhodt.sum()*dt/rhovec.sum());
-
             auto conditions = get_criticality_conditions(model, T, rhovec);
             auto z0 = rhovec[0] / rhovec.sum();
             if (z0 < 0 || z0 > 1) {
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 408b70fd0fbf277a166799ec911dd05340a0e9ad..0750c84900c7e70283dc8dc57e8e62cbcc1f12b3 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -34,6 +34,7 @@ void init_teqp(py::module& m) {
         .def_readwrite("abs_err", &TCABOptions::abs_err)
         .def_readwrite("rel_err", &TCABOptions::rel_err)
         .def_readwrite("init_dt", &TCABOptions::init_dt)
+        .def_readwrite("max_dt", &TCABOptions::max_dt)
         .def_readwrite("integration_order", &TCABOptions::integration_order)
         ;
 
diff --git a/src/tests/catch_test_multifluid.cxx b/src/tests/catch_test_multifluid.cxx
index f3180ae3dd44c7d3634c7c0b930c99c2db9a0b57..9113d36436964c2b394b403adae921c3fe8805d1 100644
--- a/src/tests/catch_test_multifluid.cxx
+++ b/src/tests/catch_test_multifluid.cxx
@@ -13,6 +13,19 @@ TEST_CASE("Trace critical locus for nitrogen + ethane", "[crit],[multifluid]")
     std::string root = "../mycp";
     const auto model = build_multifluid_model({ "Nitrogen", "Ethane" }, root, root+"/dev/mixtures/mixture_binary_pairs.json");
 
+    for (auto ifluid = 0; ifluid < 2; ++ifluid) {
+        double T0 = model.redfunc.Tc[ifluid];
+        Eigen::ArrayXd rhovec0(2); rhovec0 = 0.0; rhovec0[ifluid] = 1.0 / model.redfunc.vc[ifluid];
+
+        auto tic0 = std::chrono::steady_clock::now();
+        std::string filename = "h";
+        using ct = CriticalTracing<decltype(model), double, Eigen::ArrayXd>;
+        TCABOptions opt; opt.init_dt = 100; opt.integration_order = 1;
+        auto j = ct::trace_critical_arclength_binary(model, T0, rhovec0, filename, opt);
+        CHECK(j.size() > 3);
+        auto tic1 = std::chrono::steady_clock::now();
+    }
+    
     for (auto ifluid = 0; ifluid < 2; ++ifluid) {
         double T0 = model.redfunc.Tc[ifluid];
         Eigen::ArrayXd rhovec0(2); rhovec0 = 0.0; rhovec0[ifluid] = 1.0 / model.redfunc.vc[ifluid];