diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b6c90db4e1f5c60edf252ccf5f2e881890e75c6
--- /dev/null
+++ b/include/teqp/algorithms/VLE.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include "teqp/derivs.hpp"
+#include <Eigen/Dense>
+
+template<typename Model, typename TYPE = double>
+class IsothermPureVLEResiduals  {
+    typedef Eigen::Array<TYPE, 2, 1> EigenArray;
+    typedef Eigen::Array<TYPE, 1, 1> EigenArray1;
+    typedef Eigen::Array<TYPE, 2, 2> EigenMatrix;
+private:
+    const Model& m_model;
+    TYPE m_T;
+    EigenMatrix J;
+    EigenArray y;
+public:
+    std::size_t icall = 0;
+
+    IsothermPureVLEResiduals(const Model& model, TYPE T) : m_model(model), m_T(T) {};
+
+    const auto& get_errors() { return y; };
+
+    auto call(const EigenArray& rhovec) {
+        assert(rhovec.size() == 2);
+
+        const EigenArray1 rhovecL = rhovec.head(1);
+        const EigenArray1 rhovecV = rhovec.tail(1);
+        const auto rhomolarL = rhovecL.sum(), rhomolarV = rhovecV.sum();
+        const auto molefracs = (EigenArray1() << 1.0).finished();
+
+        using id = IsochoricDerivatives<Model,TYPE,EigenArray1>;
+        using tdx = TDXDerivatives<Model,TYPE,EigenArray1>;
+
+        const TYPE &T = m_T;
+        const TYPE R = m_model.R;
+        
+        auto derL = tdx::get_Ar0n<2>(m_model, T, rhomolarL, molefracs);
+        auto pL = rhomolarL*R*T*(1+derL[1]);
+        auto dpLdrho = R*T*(1 + 2*derL[1] + derL[2]);
+        auto hatmurL = id::build_Psir_gradient_autodiff(m_model, T, rhovecL)[0] + R*T*log(rhomolarL);
+        auto dhatmurLdrho = id::build_Psir_Hessian_autodiff(m_model, T, rhovecL)(0,0) + R*T/rhomolarL;
+
+        auto derV = tdx::get_Ar0n<2>(m_model, T, rhomolarV, molefracs);
+        auto pV = rhomolarV*R*T*(1 + derV[1]);
+        auto dpVdrho = R*T*(1 + 2*derV[1] + derV[2]);
+        auto hatmurV = id::build_Psir_gradient_autodiff(m_model, T, rhovecV)[0] + R*T*log(rhomolarV);
+        auto dhatmurVdrho = id::build_Psir_Hessian_autodiff(m_model, T, rhovecV)(0,0) + R*T/rhomolarV;
+
+        y(0) = pL - pV;
+        J(0, 0) = dpLdrho;
+        J(0, 1) = -dpVdrho;
+
+        y(1) = hatmurL - hatmurV;
+        J(1, 0) = dhatmurLdrho;
+        J(1, 1) = -dhatmurVdrho;
+
+        icall++;
+        return y;
+    }
+    auto Jacobian(const EigenArray& rhovec){
+        return J;
+    }
+};
diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx
index 52e9b90e90f2365b7ef264bf3559f534f6ea0d8d..7e38397af8fb4fea615d52dc7e4bed7f864c752b 100644
--- a/src/tests/catch_tests.cxx
+++ b/src/tests/catch_tests.cxx
@@ -3,6 +3,7 @@
 
 #include "teqp/core.hpp"
 #include "teqp/models/pcsaft.hpp"
+#include "teqp/algorithms/VLE.hpp"
 
 auto build_vdW_argon() {
     double Omega_b = 1.0 / 8, Omega_a = 27.0 / 64;
@@ -228,4 +229,17 @@ TEST_CASE("Test psir gradient", "") {
     CAPTURE(err1);
     CHECK(err0 < 1e-12);
     CHECK(err1 < 1e-12);
+}
+
+TEST_CASE("Test pure VLE", "") {
+    const auto model = build_vdW_argon();
+    double T = 100.0;
+    auto resid = IsothermPureVLEResiduals(model, T);
+    auto rhovec = (Eigen::ArrayXd(2) << 22834.056386882046, 1025.106554560764).finished();
+    auto r0 = resid.call(rhovec);
+    auto J = resid.Jacobian(rhovec);
+    auto v = J.matrix().colPivHouseholderQr().solve(-r0.matrix()).array().eval();
+    auto rhovec1 = rhovec + v;
+    auto r1 = resid.call(rhovec1);
+    CHECK((r0.cwiseAbs() > r1.cwiseAbs()).eval().all());
 }
\ No newline at end of file