diff --git a/externals/mcx b/externals/mcx
index 2a2e79661b07e884793b2e0831926692d313b1ff..414a4a6677b3923064c8ea20ce84f11cbc0b7b10 160000
--- a/externals/mcx
+++ b/externals/mcx
@@ -1 +1 @@
-Subproject commit 2a2e79661b07e884793b2e0831926692d313b1ff
+Subproject commit 414a4a6677b3923064c8ea20ce84f11cbc0b7b10
diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index eb71882d064e7871ca6cde51178b38e9aa0835c0..efb110e48feca58bb6726927190094d006c1cd97 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -132,6 +132,7 @@ auto get_derivs(const Model& model, const TType T, const RhoType& rhovec) {
     auto derivs = derivatives(wrapper, wrt(varsigma, varsigma, varsigma, varsigma), at(varsigma));
     auto psir_derivs = std::valarray<double>(&derivs[0], derivs.size());
 #else
+    using namespace mcx;
     // Calculate the first through fourth derivative of Psi^r w.r.t. sigma_1
     std::valarray<MultiComplex<double>> v0(ei.v0.size()); for (auto i = 0; i < ei.v0.size(); ++i) { v0[i] = ei.v0[i]; }
     std::valarray<MultiComplex<double>> rhovecmcx(rhovec.size());  for (auto i = 0; i < rhovec.size(); ++i) { rhovecmcx[i] = rhovec[i]; }
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 0ad082427e2ca948468af1751246b87a763e3803..51dd9506df418875efcc98efb9a239dadcbe07c4 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -33,8 +33,8 @@ typename ContainerType::value_type derivT(const FuncType& f, TType T, const Cont
 */
 template <typename TType, typename ContainerType, typename FuncType>
 typename ContainerType::value_type derivTmcx(const FuncType& f, TType T, const ContainerType& rho) {
-    using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>;
-    fcn_t wrapper = [&rho, &f](const MultiComplex<TType>& T_) {return f(T_, rho); };
+    using fcn_t = std::function<mcx::MultiComplex<double>(const mcx::MultiComplex<double>&)>;
+    fcn_t wrapper = [&rho, &f](const mcx::MultiComplex<TType>& T_) {return f(T_, rho); };
     auto ders = diff_mcx1(wrapper, T, 1);
     return ders[0];
 }
@@ -113,6 +113,7 @@ typename ContainerType::value_type get_B2vir(const Model& model, const TType T,
 template <typename Model, typename TType, typename ContainerType>
 auto get_Bnvir(const Model& model, int Nderiv, const TType T, const ContainerType& molefrac) {
     // B_n = lim_rho\to 0 d^{n-1}alphar/drho^{n-1}|T,z
+    using namespace mcx;
     using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>;
     fcn_t f = [&model, &T, &molefrac](const MultiComplex<double>& rho_) -> MultiComplex<double> { return model.alphar(T, rho_, molefrac); };
     return diff_mcx1(f, 0.0, Nderiv, true /* and_val */);
@@ -147,16 +148,7 @@ typename ContainerType::value_type get_pr(const Model& model, const TType T, con
     return get_Ar01(model, T, rhotot_, rhovec / rhotot_) * rhotot_ * model.R * T;
 }
 
-// Generic setting functions to handle Eigen types and STL types with the same interface
-template<typename MatrixLike, typename Integer, typename ValType>
-void setval(MatrixLike &m, Integer i, Integer j, const ValType val) {
-    m(i,j) = val;
-}
 
-// Partial specialization for valarray "matrix"
-template <> void setval<std::valarray<std::valarray<double>>, std::size_t, double>(std::valarray<std::valarray<double>>& m, std::size_t i, std::size_t j, const double val) {
-    m[i][j] = val;
-}
 
 /***
 * \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
@@ -188,27 +180,14 @@ template<typename Model, typename TType, typename RhoType>
 auto build_Psir_Hessian_mcx(const Model& model, const TType &T, const RhoType& rho) {
     // Double derivatives in each component's concentration
     // N^N matrix (symmetric)
+    using namespace mcx;
 
     // Lambda function for getting Psir with multicomplex concentrations
-    auto func = [&model, &T](const std::vector<MultiComplex<double>>& rhovec) {
-        std::valarray<MultiComplex<double>> xs(&(rhovec[0]), rhovec.size());
-        return get_Psir(model, T, xs);
+    using fcn_t = std::function< MultiComplex<double>(const std::valarray<MultiComplex<double>>&)>;
+    fcn_t func = [&model, &T](const std::valarray<MultiComplex<double>>& rhovec) -> MultiComplex<double> {
+        return get_Psir(model, T, rhovec);
     };
-    // The set of values around which the pertubations will happen
-    const std::size_t N = rho.size();
-    std::vector<double> xs(std::begin(rho), std::end(rho));
-
-    Eigen::MatrixXd H(N, N);
-    
-    for (std::size_t i = 0; i < rho.size(); ++i) {
-        for (std::size_t j = i; j < rho.size(); ++j) {
-            std::vector<int> order = { 0, 0 };
-            order[i] += 1;
-            order[j] += 1;
-            auto val = diff_mcxN<double>(func, xs, order);
-            setval(H,i,j,val);
-            setval(H,j,i,val);
-        }
-    }
+    using mattype = Eigen::ArrayXXd;
+    auto H = get_Hessian<mattype, fcn_t, std::valarray<double>, HessianMethods::Multiple>(func, rho);
     return H;
 }
\ No newline at end of file