diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 81276a0..65c46fd 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,38 +1,46 @@
 set(PROJECT_EXAMPLE_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 
 
 # ========== SpecMiCP =====================================================
 
 # ---------- adimensional system ------------------------------------------
 
 #            thermocarbo
 
 add_executable(ex_adim_thermocarbo
     ${PROJECT_EXAMPLE_DIR}/specmicp/adimensional/thermocarbo.cpp
 )
 
 target_link_libraries(ex_adim_thermocarbo specmicp specmicp_database)
 
+#            equilibrium curve
+
+add_executable(ex_adim_equilibriumcurve
+    ${PROJECT_EXAMPLE_DIR}/specmicp/adimensional/equilibrium_curve.cpp
+)
+
+target_link_libraries(ex_adim_equilibriumcurve specmicp specmicp_database)
+
 # ========== ReactMiCP ====================================================
 
 # ---------- leaching -----------------------------------------------------
 
 add_executable(leaching 
     ${PROJECT_EXAMPLE_DIR}/reactmicp/leaching.cpp
 )
 target_link_libraries(leaching reactmicp specmicp specmicp_database)
 
 # ---------- carbonation --------------------------------------------------
 
 add_executable(carbonation
     ${PROJECT_EXAMPLE_DIR}/reactmicp/carbonation.cpp
 )
 target_link_libraries(carbonation reactmicp specmicp specmicp_database)
 
 # ---------- diffusion ----------------------------------------------------
 
 add_executable(diffusion
     ${PROJECT_EXAMPLE_DIR}/reactmicp/diffusion.cpp
 )
 target_link_libraries(diffusion reactmicp specmicp specmicp_database)
 
diff --git a/examples/specmicp/adimensional/equilibrium_curve.cpp b/examples/specmicp/adimensional/equilibrium_curve.cpp
new file mode 100644
index 0000000..e714461
--- /dev/null
+++ b/examples/specmicp/adimensional/equilibrium_curve.cpp
@@ -0,0 +1,139 @@
+#include <iostream>
+
+#include "utils/log.hpp"
+#include "specmicp/adimensional/adimensional_system_solver.hpp"
+#include "specmicp/adimensional/adimensional_system_solution.hpp"
+#include "specmicp/problem_solver/formulation.hpp"
+#include "specmicp/problem_solver/dissolver.hpp"
+#include "specmicp/adimensional/adimensional_system_solution_extractor.hpp"
+
+#include "specmicp/adimensional/equilibrium_curve.hpp"
+
+class LeachingEquilibriumCurve: public specmicp::EquilibriumCurve
+{
+public:
+    specmicp::scalar_t step = 25;
+
+    LeachingEquilibriumCurve()
+    {
+        specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
+        std::map<std::string, std::string> swapping ({
+                                                  {"H[+]","HO[-]"},
+                                                  {"Si(OH)4", "SiO(OH)3[-]"},
+                                                    });
+        thedatabase.swap_components(swapping);
+        thedatabase.remove_gas_phases();
+        specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
+        set_database(raw_data);
+
+        specmicp::Formulation formulation;
+        specmicp::scalar_t mult = 7e3;
+
+        specmicp::scalar_t m_c3s = mult*0.7;
+        specmicp::scalar_t m_c2s = mult*0.3;
+        specmicp::scalar_t wc = 0.5;
+        specmicp::scalar_t m_water = wc*1e-3*(
+                      m_c3s*(3*56.08+60.08)
+                    + m_c2s*(2*56.06+60.08)
+                    );
+
+        formulation.mass_solution = m_water;
+        formulation.amount_minerals = {
+            {"C3S", m_c3s},
+            {"C2S", m_c2s},
+        };
+
+        specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
+        id_h2o = thedatabase.component_label_to_id("H2O");
+        id_ho  = thedatabase.component_label_to_id("HO[-]");
+        id_ca = thedatabase.component_label_to_id("Ca[2+]");
+
+        conditions() = specmicp::AdimensionalSystemBC(total_concentrations);
+        conditions().charge_keeper = id_ho;
+
+        specmicp::AdimensionalSystemSolverOptions& options = solver_options();
+        options.solver_options.maxstep = 10.0;
+        options.solver_options.max_iter = 100;
+        options.solver_options.maxiter_maxstep = 100;
+        options.solver_options.use_crashing = false;
+        options.solver_options.use_scaling = false;
+        options.solver_options.factor_descent_condition = -1;
+        options.solver_options.factor_gradient_search_direction = 100;
+        options.solver_options.projection_min_variable = 1e-9;
+        options.solver_options.fvectol = 1e-6;
+        options.solver_options.steptol = 1e-14;
+        options.system_options.non_ideality_tolerance = 1e-10;
+
+        solve_first_problem();
+
+        std::cout << "total concentration \t TOTaq \t TOTs \t S^t_m \t pH";
+        for (specmicp::index_t mineral: raw_data->range_mineral())
+        {
+            std::cout << "\t" << raw_data->labels_minerals[mineral];
+        }
+        std::cout << std::endl;
+
+
+        specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set);
+        specmicp::scalar_t tot_conc_ca = conditions().total_concentrations(id_ca);
+        total_steps = std::ceil(tot_conc_ca / step) + 1;
+        cnt = 0;
+        solution = specmicp::Matrix::Zero(total_steps, 2);
+
+
+        output();
+
+    }
+
+    void output()
+    {
+        specmicp::AdimensionalSystemSolutionExtractor sol(current_solution(), database(), solver_options().units_set);
+        std::cout << conditions().total_concentrations(id_ca) << "\t"
+                  << sol.total_aqueous_concentration(id_ca)   << "\t"
+                  << sol.total_solid_concentration(id_ca) << "\t"
+                  << sol.total_saturation_minerals() << "\t"
+                  << sol.pH();
+        for (specmicp::index_t mineral: database()->range_mineral())
+        {
+            std::cout << "\t" << sol.mole_concentration_mineral(mineral);
+        }
+        std::cout << std::endl;
+        solution(cnt, 0) = sol.total_solid_concentration(id_ca);
+        solution(cnt, 1) = sol.total_aqueous_concentration(id_ca);
+
+        ++cnt;
+    }
+
+    specmicp::index_t nb_steps() const
+    {
+        return total_steps;
+    }
+
+    void update_problem()
+    {
+        conditions().total_concentrations(id_ca) -=   step;
+        conditions().total_concentrations(id_ho) += 2*step;
+    }
+
+    specmicp::Matrix& equilibrium_curve() {return solution;}
+
+private:
+    mutable specmicp::index_t cnt;
+    specmicp::index_t total_steps;
+
+    mutable specmicp::Matrix solution;
+    specmicp::index_t id_h2o;
+    specmicp::index_t id_ho;
+    specmicp::index_t id_ca;
+};
+
+
+int main()
+{
+    LeachingEquilibriumCurve solver;
+    for (int i=1; i<solver.nb_steps(); ++i)
+    {
+        solver.run_step();
+    }
+    std::cout << solver.equilibrium_curve() << std::endl;
+}