diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 61b0f627b21cc67f5bb8af76cba983181a85ff56..6726bab6925af0d4e155db4cbf3678bd88ac8418 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -6,8 +6,10 @@ caller(const FuncType& f, TType T, const ContainerType& rho) {
     return f(T, rho);
 }
 
-/// Given a function, use complex step derivatives to calculate the derivative with respect to the first variable
-/// which here is temperature
+/***
+* \brief Given a function, use complex step derivatives to calculate the derivative with 
+* respect to the first variable which here is temperature
+*/
 template <typename TType, typename ContainerType, typename FuncType>
 typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
 derivT(const FuncType& f, TType T, const ContainerType& rho) {
@@ -15,7 +17,9 @@ derivT(const FuncType& f, TType T, const ContainerType& rho) {
     return f(std::complex<TType>(T, h), rho).imag() / h;
 }
 
-/// Given a function, use complex step derivatives to calculate the derivative with respect to the given composition variable
+/***
+* \brief Given a function, use complex step derivatives to calculate the derivative with respect to the given composition variable
+*/
 template <typename TType, typename ContainerType, typename FuncType, typename Integer>
 typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
 derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer i) {
@@ -29,6 +33,9 @@ derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer i) {
     return f(T, rhocom).imag() / h;
 }
 
+/***
+* \brief Calculate the Psir=ar*rho
+*/
 template <typename TType, typename ContainerType, typename Model>
 typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
 get_Psir(const Model& model, const TType T, const ContainerType& rhovec) {
@@ -37,12 +44,15 @@ get_Psir(const Model& model, const TType T, const ContainerType& rhovec) {
     return model.alphar(T, rhovec)*model.R*T*rhotot_;
 }
 
-/**
-/// Calculate the residual pressure from derivatives of alphar
+/***
+* \brief Calculate the residual pressure from derivatives of alphar
 */
 template <typename Model, typename TType, typename ContainerType>
-typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
-get_pr(const Model& model, const TType T, const ContainerType& rhovec) {
+typename ContainerType::value_type get_pr(
+    const Model& model, 
+    const TType T, 
+    const ContainerType& rhovec) 
+{
     using container = decltype(rhovec);
     auto rhotot_ = std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
     decltype(rhovec[0] * T) pr = 0.0;
@@ -52,17 +62,33 @@ get_pr(const Model& model, const TType T, const ContainerType& rhovec) {
     return pr*rhotot_*model.R*T;
 }
 
-/**
-/// Calculate the residual entropy (s^+=-sr/R) from derivatives of alphar
+/***
+* \brief Calculate the residual entropy (s^+ = -sr/R) from derivatives of alphar
 */
 template <typename Model, typename TType, typename ContainerType>
-typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
-get_splus(const Model& model, const TType T, const ContainerType& rhovec) {
+typename ContainerType::value_type get_splus(
+    const Model& model, 
+    const TType T, 
+    const ContainerType& rhovec)
+{
     return model.alphar(T, rhovec) + T*derivT([&model](const auto& T, const auto& rhovec) { return model.alphar(T, rhovec); }, T, rhovec);
 }
 
+/***
+* \brief Calculate the Hessian of Psir = ar*rho w.r.t. the molar concentrations
+* 
+* Requires the use of multicomplex derivatives to calculate second partial derivatives
+*/
 template<typename Model, typename TType, typename RhoType>
 auto build_Psir_Hessian(const Model& model, const TType T, const RhoType& rho) {
     // Double derivatives in each component's concentration
-
+    // N^N matrix (symmetric)
+    for (auto i = 0; i < rho.size(); ++i) {
+        for (auto j = i; j < rho.size(); ++j) {
+            auto val = diff_mcxN();
+            H(i,j) = val;
+            H(j,i) = val;
+        }
+    }
+    return H;
 }
\ No newline at end of file