diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index bec942c..53a0c4b 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,53 +1,53 @@
 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)
 
 # ---------- EquilibriumCurve ----------------------------------------------------
 
 add_executable(ex_reactmicp_equilibriumcurve
     ${PROJECT_EXAMPLE_DIR}/reactmicp/equilibrium_curve.cpp
 )
-target_link_libraries(ex_reactmicp_equilibriumcurve reactmicp specmicp specmicp_database)
+target_link_libraries(ex_reactmicp_equilibriumcurve reactmicp dfpm specmicp specmicp_database)
 
diff --git a/examples/reactmicp/equilibrium_curve.cpp b/examples/reactmicp/equilibrium_curve.cpp
index 0c49a7d..8e12fe1 100644
--- a/examples/reactmicp/equilibrium_curve.cpp
+++ b/examples/reactmicp/equilibrium_curve.cpp
@@ -1,148 +1,367 @@
 #include <iostream>
 #include "utils/log.hpp"
 #include "reactmicp/equilibrium_curve/chemistry.hpp"
 #include "specmicp/problem_solver/dissolver.hpp"
 #include "specmicp/problem_solver/formulation.hpp"
 
+#include "reactmicp/equilibrium_curve/eqcurve_extractor.hpp"
+#include "reactmicp/equilibrium_curve/eqcurve_coupler.hpp"
+#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
+#include "dfpm/meshes/uniform_mesh1d.hpp"
 
-void test_chemistry() {
+#include "reactmicp/equilibrium_curve/eqcurve_solid_transport.hpp"
+
+specmicp::Matrix test_chemistry() {
 
     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();
 
     specmicp::Formulation formulation;
-    specmicp::scalar_t mult = 7e3;
+    specmicp::scalar_t mult = 6.5e3;
 
     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);
     specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
     specmicp::index_t id_ho  = thedatabase.component_label_to_id("HO[-]");
     specmicp::index_t id_ca = thedatabase.component_label_to_id("Ca[2+]");
 
     specmicp::AdimensionalSystemBC conditions(total_concentrations);
     conditions.charge_keeper = id_ho;
 
     specmicp::AdimensionalSystemSolverOptions 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;
 
     specmicp::reactmicp::eqcurve::EquilibriumCurveSpeciation spec_solver(raw_data, conditions, id_ca, options);
 
-    std::cout << spec_solver.get_equilibrium_curve(500, 1.0) << std::endl;
+    return spec_solver.get_equilibrium_curve(0.05, -500);
 }
 
 
-void test_chemistry_with_al()
+specmicp::Matrix test_chemistry_with_al()
 {
 
     specmicp::database::Database thedatabase("../data/cemdata_specmicp.js");
     std::map<std::string, std::string> swapping ({
                                               {"H[+]","HO[-]"},
                                               {"Si(OH)4", "SiO(OH)3[-]"},
                                               {"Al[3+]","Al(OH)4[-]"}
                                                 });
     thedatabase.swap_components(swapping);
     thedatabase.remove_gas_phases();
     specmicp::RawDatabasePtr raw_data = thedatabase.get_database();
 
     specmicp::Formulation formulation;
-    specmicp::scalar_t mult = 6e3;
+    specmicp::scalar_t mult = 6.5e3;
 
     specmicp::scalar_t m_c3s = mult*0.6;
     specmicp::scalar_t m_c2s = mult*0.2;
     specmicp::scalar_t m_c3a = mult*0.10;
     specmicp::scalar_t m_gypsum = mult*0.10;
     specmicp::scalar_t wc = 0.8;
     specmicp::scalar_t m_water = wc*1e-3*(
                   m_c3s*(3*56.08+60.08)
                 + m_c2s*(2*56.06+60.08)
                 + m_c3a*(3*56.08+101.96)
                 + m_gypsum*(56.08+80.06+2*18.02)
                 );
 
     formulation.mass_solution = m_water;
     formulation.amount_minerals = {
         {"C3S", m_c3s},
         {"C2S", m_c2s},
         {"C3A", m_c3a},
         {"Gypsum", m_gypsum}
     };
     formulation.minerals_to_keep = {
         "Portlandite",
         "CSH,jennite",
         "CSH,tobermorite",
         "SiO2,am",
         "Al(OH)3,am",
         "Monosulfoaluminate",
         "Straetlingite",
         "Gypsum",
         "Ettringite",
     };
 
     for (specmicp::index_t component: raw_data->range_component())
     {
         std::cout << raw_data->labels_basis[component] << std::endl;
     }
 
     specmicp::Vector total_concentrations = specmicp::Dissolver(raw_data).dissolve(formulation);
-    specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
+    //specmicp::index_t id_h2o = thedatabase.component_label_to_id("H2O");
     specmicp::index_t id_ho  = thedatabase.component_label_to_id("HO[-]");
     specmicp::index_t id_ca = thedatabase.component_label_to_id("Ca[2+]");
 
     specmicp::AdimensionalSystemBC conditions(total_concentrations);
     conditions.charge_keeper = id_ho;
 
     specmicp::AdimensionalSystemSolverOptions options;
     options.solver_options.maxstep = 20.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;
 
     specmicp::reactmicp::eqcurve::EquilibriumCurveSpeciation spec_solver(raw_data, conditions, id_ca, options);
 
-    std::cout << spec_solver.get_equilibrium_curve(10.0, -750.0) << std::endl;
+    return spec_solver.get_equilibrium_curve(0.05, -500.0);
+}
+
+void test_eqcurve_extractor()
+{
+    specmicp::reactmicp::eqcurve::EquilibriumCurveExtractor extract(test_chemistry_with_al());
+
+
+    specmicp::index_t index = extract.find_point(111.0);
+    std::cout << "111.0 \t " << extract.totsolid_concentration(index) << "\t"
+              << extract.totaq_concentration(index) << "\t"
+              << extract.porosity(index) << "\t"
+              << extract.diffusion_coefficient(index) << std::endl;
+}
+
+void test_diffeqcurve()
+{
+    specmicp::Matrix eqcurve = test_chemistry_with_al();
+    eqcurve.col(0) *= 1e-6; //mol/m3 -> mol/cm3
+    eqcurve.col(1) *= 1e-3; //mol/kg -> mol/cm3
+
+    std::cout << eqcurve << std::endl;
+
+    specmicp::scalar_t radius = 3.5; //cm
+    specmicp::scalar_t height = 8.0; //cm
+
+    specmicp::scalar_t dx = 0.0025;
+    specmicp::index_t additional_nodes = 1;
+
+    radius = radius + additional_nodes*dx;
+    specmicp::index_t nb_nodes =10+additional_nodes;
+
+    specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
+
+
+    specmicp::dfpmsolver::ParabolicDriverOptions options;
+    options.step_tolerance = 1e-10;
+    options.residuals_tolerance = 1e-8;
+    options.quasi_newton = 1;
+
+
+    specmicp::reactmicp::eqcurve::EquilibriumCurveCoupler solver(eqcurve, the_mesh, options);
+
+
+    specmicp::scalar_t sum_0 =0;
+    for (specmicp::index_t node=0; node< the_mesh->nb_nodes(); ++node)
+    {
+        sum_0 += solver.solid_concentrations()(node)*the_mesh->get_volume_cell(node);
+        std::cout << the_mesh->get_volume_cell(node) << std::endl;
+    }
+    specmicp::scalar_t dt = 0.5;
+    specmicp::scalar_t total = 0;
+
+    std::cout << total << "\t" << 0.0 << "\t" << sum_0 << "\t" << 0.0 << std::endl;
+    specmicp::index_t k=0;
+    while (total < 0.4)
+    {
+        solver.run_step(dt);
+        total += dt/3600/24;
+        if (k % 5000 == 0)
+        {
+            specmicp::scalar_t sum =0;
+            for (specmicp::index_t node=0; node< the_mesh->nb_nodes(); ++node)
+            {
+                sum += solver.solid_concentrations()(node)*the_mesh->get_volume_cell(node);
+            }
+            std::cout << total << "\t" << std::sqrt(total) << "\t" << sum << "\t" << (sum_0 -sum)/(1.75929e-2) << std::endl;
+
+        }
+        ++k;
+    }
+    //std::cout << solver.solid_concentrations() << std::endl;
+}
+
+void test_eqcurve_solid()
+{
+    //specmicp::Matrix eq_curve = test_chemistry();
+    specmicp::Matrix eq_curve = test_chemistry_with_al();
+    eq_curve.col(0) *= 1e-6; //mol/m3 -> mol/cm3
+    eq_curve.col(1) *= 1e-3; //mol/kg -> mol/cm3
+    for (specmicp::index_t ind=1; ind<eq_curve.rows(); ++ind)
+    {
+        if (eq_curve(ind, 1) >= eq_curve(ind-1, 1))
+            eq_curve(ind, 1) = eq_curve(ind-1, 1);
+    }
+
+
+    std::cout << eq_curve << std::endl;
+
+    specmicp::scalar_t radius = 3.5; //cm
+    specmicp::scalar_t height = 8.0; //cm
+
+    specmicp::scalar_t dx = 0.005;
+    specmicp::index_t additional_nodes = 1;
+
+    radius = radius + additional_nodes*dx;
+    specmicp::index_t nb_nodes = 50+additional_nodes;
+
+    specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::axisymmetric_uniform_mesh1d(nb_nodes, radius, dx, height);
+    //specmicp::mesh::Mesh1DPtr the_mesh = specmicp::mesh::uniform_mesh1d(nb_nodes, dx, 5);
+
+    specmicp::dfpmsolver::ParabolicDriverOptions options;
+    options.step_tolerance = 1e-12;
+    options.residuals_tolerance = 1e-6;
+    options.sparse_solver = specmicp::SparseSolver::GMRES;
+    //options.linesearch = specmicp::dfpmsolver::ParabolicLinesearch::Strang;
+    options.alpha = 1.0;
+    options.quasi_newton = 1;
+    options.maximum_step_length = 10;
+
+
+    specmicp::reactmicp::eqcurve::SolidDiffusion program(the_mesh, eq_curve, {0,});
+    specmicp::dfpmsolver::ParabolicDriver< specmicp::reactmicp::eqcurve::SolidDiffusion> solver(program);
+    solver.get_options() = options;
+
+    solver.set_scaling(specmicp::Vector::Constant(program.get_neq(), 1e6));
+
+    specmicp::Vector variables(nb_nodes);
+    variables(0) = eq_curve(eq_curve.rows()-10, 0);
+    for (specmicp::index_t node=1; node<the_mesh->nb_nodes(); ++node)
+    {
+        variables(node) = eq_curve(10, 0);
+    }
+
+
+    specmicp::scalar_t sum_0 =0;
+    for (specmicp::index_t node=0; node< the_mesh->nb_nodes(); ++node)
+    {
+        sum_0 += variables(node)*the_mesh->get_volume_cell(node);
+        std::cout << the_mesh->get_volume_cell(node) << std::endl;
+    }
+    specmicp::scalar_t dt = 10.0;
+    specmicp::scalar_t total = 0;
+
+    std::cout << total << "\t" << 0.0 << "\t" << sum_0 << "\t" << 0.0 << std::endl;
+    specmicp::index_t k=0;
+    while (total < 50)
+    {
+        //std::cout << " ==== TIMESTEP === " << std::endl;
+        solver.solve_timestep(dt, variables);
+        for (int node=0; node<the_mesh->nb_nodes(); ++node)
+           if (variables(node) < 1e-6) variables(node) = 0;
+        //std::cout << solver.get_perfs().nb_iterations << std::endl;
+        total += dt/3600/24;
+        if (k % 5000 == 0)
+        {
+            specmicp::scalar_t sum =0;
+            for (specmicp::index_t node=0; node< the_mesh->nb_nodes(); ++node)
+            {
+                sum += variables(node)*the_mesh->get_volume_cell(node);
+            }
+            std::cout << total << "\t" << std::sqrt(total) << "\t" << sum << "\t" << (sum_0 -sum)/(1.75929e-2) << std::endl;
+
+        }
+        ++k;
+    }
+    std::cout << variables << std::endl;
+}
+
+void test_interpolator()
+{
+    specmicp::Matrix mat(5,4);
+
+    mat << 1, 1, 1, 1,
+           2, 1, 2, 0,
+           3, 1, 3, -1,
+           4, 1, 4, -2,
+           5, 1, 5, -3;
+
+    specmicp::reactmicp::eqcurve::EquilibriumCurveExtractor interpolator(mat);
+
+    std::cout << " " << interpolator.slope(0, 1) << " ?== " << 0 << std::endl;
+    std::cout << " " << interpolator.slope(0, 2) << " ?== " << 1 << std::endl;
+    std::cout <<        interpolator.slope(0, 3) << " ?== " << -1 << std::endl;
+
+    std::cout << interpolator.find_point(1.5)  << " ? == " << 0 << std::endl;
+    std::cout << interpolator.find_point(3.5)  << " ? == " << 2 << std::endl;
+
+    std::cout << interpolator.interpolate(2, 3.5, 2) << " ? == " <<  3.5 << std::endl;
+    std::cout << interpolator.interpolate(2, 3.5, 3) << " ? == " << -1.5 << std::endl;
+
+    specmicp::Matrix mat2(5,4);
+
+    mat2 << 5, 1, 1, 1,
+           4, 1, 2, 0,
+           3, 1, 3, -1,
+           2, 1, 4, -2,
+           1, 1, 5, -3;
+
+    specmicp::reactmicp::eqcurve::EquilibriumCurveExtractor interpolator2(mat2);
+
+    std::cout << " " << interpolator2.slope(0, 1) << " ?== " <<  0 << std::endl;
+    std::cout << " " << interpolator2.slope(0, 2) << " ?== " << -1 << std::endl;
+    std::cout <<        interpolator2.slope(0, 3) << " ?== " << +1 << std::endl;
+
+    std::cout << interpolator2.find_point(1.5)  << " ? == " << 3 << std::endl;
+    std::cout << interpolator2.find_point(3.5)  << " ? == " << 1 << std::endl;
+
+    std::cout << interpolator2.interpolate(1, 3.5, 2) << " ? == " <<  2.5 << std::endl;
+    std::cout << interpolator2.interpolate(1, 3.5, 3) << " ? == " << -0.5 << std::endl;
+
+    std::cout << interpolator2.interpolate(4, 1.0, 2) << " ? == " <<  5 << std::endl;
+    std::cout << interpolator2.interpolate(4, 1.0, 3) << " ? == " << -3 << std::endl;
+
 }
 
 int main()
 {
     specmicp::stdlog::ReportLevel() = specmicp::logger::Warning;
     specmicp::logger::ErrFile::stream() = &std::cerr;
 
     //test_chemistry();
-    test_chemistry_with_al();
+    //std::cout <<  test_chemistry_with_al() << std::endl;
+
+    //test_eqcurve_extractor();
+
+    //test_diffeqcurve();
+
+    test_eqcurve_solid();
+
+    test_interpolator();
+
+    return EXIT_SUCCESS;
 }