diff --git a/src/BaseSolver.cpp b/src/BaseSolver.cpp
index 05c837e..2ba46cb 100644
--- a/src/BaseSolver.cpp
+++ b/src/BaseSolver.cpp
@@ -140,12 +140,11 @@ Eigen::VectorXi BaseSolver::extract_slack_bus_id(const Eigen::VectorXi & pv,
     // pv: list of index of pv nodes
     // pq: list of index of pq nodes
     // nb_bus: total number of bus in the grid
-    // returns: res: the id of the unique slack bus (throw an error if no slack bus is found)
-    // /!\ does not support multiple slack bus!!!
+    // returns: res: the ids of all the slack buses (by def: not PV and not PQ)
 
     Eigen::VectorXi res(nb_bus - pv.size() - pq.size());
     Eigen::Index i_res = 0;
-    bool found=false;
+    
     // run through both pv and pq nodes and declare they are not slack bus
     std::vector<bool> tmp(nb_bus, true);
     for(unsigned int k=0; k < pv.size(); ++k)
@@ -163,16 +162,10 @@ Eigen::VectorXi BaseSolver::extract_slack_bus_id(const Eigen::VectorXi & pv,
         {
             res[i_res] = k;
             ++i_res;
-            // if(!found){
-            //     res = k;
-            //     found = true;
-            // }else{
-            //     throw std::runtime_error("BaseSolver::extract_slack_bus_id: multiple slack bus found on your grid !");
-            // }
         }
     }
-    // if(res == -1){
     if(res.size() != i_res){
+        // TODO DEBUG MODE
         throw std::runtime_error("BaseSolver::extract_slack_bus_id: No slack bus is found in your grid");
     }
     return res;
diff --git a/src/CKTSOSolver.cpp b/src/CKTSOSolver.cpp
index 0b673f7..8d5b734 100644
--- a/src/CKTSOSolver.cpp
+++ b/src/CKTSOSolver.cpp
@@ -11,6 +11,8 @@
 #include "CKTSOSolver.h"
 #include <iostream>
 
+const bool CKTSOLinearSolver::CAN_SOLVE_MAT = false;
+
 
 ErrorType CKTSOLinearSolver::reset(){
     // free everything
diff --git a/src/CKTSOSolver.h b/src/CKTSOSolver.h
index 101cb62..2ed7e34 100644
--- a/src/CKTSOSolver.h
+++ b/src/CKTSOSolver.h
@@ -64,6 +64,9 @@ class CKTSOLinearSolver
         ErrorType initialize(Eigen::SparseMatrix<real_type> & J);
         ErrorType solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized);
 
+        // can this linear solver solve problem where RHS is a matrix
+        static const bool CAN_SOLVE_MAT;
+        
         // prevent copy and assignment
         CKTSOLinearSolver(const CKTSOLinearSolver & other) = delete;
         CKTSOLinearSolver & operator=( const CKTSOLinearSolver & ) = delete;
diff --git a/src/DCSolver.h b/src/DCSolver.h
index e8ee780..7b20b04 100644
--- a/src/DCSolver.h
+++ b/src/DCSolver.h
@@ -43,11 +43,25 @@ class BaseDCSolver: public BaseSolver
         BaseDCSolver( const BaseSolver & ) =delete ;
         BaseDCSolver & operator=( const BaseSolver & ) =delete;
 
+    protected:
+        void fill_mat_bus_id(int nb_bus_solver);
+        void fill_dcYbus_noslack(int nb_bus_solver, const Eigen::SparseMatrix<cplx_type> & ref_mat);
+
+        // remove_slack_buses: res_mat is initialized and make_compressed in this function
+        template<typename ref_mat_type>  // ref_mat_type should be `real_type` or `cplx_type`
+        void remove_slack_buses(int nb_bus_solver, const Eigen::SparseMatrix<ref_mat_type> & ref_mat, Eigen::SparseMatrix<real_type> & res_mat);
+
     protected:
         LinearSolver  _linear_solver;
         bool need_factorize_;
-        RealVect Sbus_noslack_;
-        Eigen::SparseMatrix<real_type> Ybus_noslack_;
+
+        // save this not to recompute them when not needed
+        RealVect dcSbus_noslack_;
+        Eigen::SparseMatrix<real_type> dcYbus_noslack_;
+        Eigen::VectorXi my_pv_;
+        Eigen::VectorXi slack_buses_ids_solver_;
+        // -1 if bus is slack , else the id of the row / column used in the linear solver representing this bus
+        Eigen::VectorXi mat_bus_id_;   // formerly `ybus_to_me`
 
 };
 
diff --git a/src/DCSolver.tpp b/src/DCSolver.tpp
index 114f76c..e0b1ad1 100644
--- a/src/DCSolver.tpp
+++ b/src/DCSolver.tpp
@@ -33,47 +33,24 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     #ifdef __COUT_TIMES
         auto timer_preproc = CustTimer();
     #endif // __COUT_TIMES
-    Eigen::SparseMatrix<real_type> dcYbus = Eigen::SparseMatrix<real_type>(nb_bus_solver - 1, nb_bus_solver - 1);
 
     const CplxVect & Sbus_tmp = Sbus;
 
     // TODO SLACK (for now i put all slacks as PV, except the first one)
     // this should be handled in Sbus, because we know the amount of power absorbed by the slack
     // so we can compute it correctly !
-    const Eigen::VectorXi my_pv = retrieve_pv_with_slack(slack_ids, pv);
+    my_pv_ = retrieve_pv_with_slack(slack_ids, pv);
     // const Eigen::VectorXi & my_pv = pv;
 
     // find the slack buses
-    Eigen::VectorXi slack_bus_ids_solver = extract_slack_bus_id(my_pv, pq, nb_bus_solver);
+    slack_buses_ids_solver_ = extract_slack_bus_id(my_pv_, pq, nb_bus_solver);
+
     // corresp bus -> solverbus
-    Eigen::VectorXi ybus_to_me = Eigen::VectorXi::Constant(nb_bus_solver, -1);
-    // Eigen::VectorXi me_to_ybus = Eigen::VectorXi::Constant(nb_bus_solver - slack_bus_ids_solver.size(), -1);
-    int solver_id = 0;
-    for (int ybus_id=0; ybus_id < nb_bus_solver; ++ybus_id){
-        if(isin(ybus_id, slack_bus_ids_solver)) continue;  // I don't add anything to the slack bus
-        ybus_to_me(ybus_id) = solver_id;
-        // me_to_ybus(solver_id) = ybus_id;
-        ++solver_id;
-    }
+    fill_mat_bus_id(nb_bus_solver);
+
     // remove the slack bus from Ybus
     // and extract only real part
-    // TODO see if "prune" might work here https://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html#title29
-    std::vector<Eigen::Triplet<real_type> > tripletList;
-    tripletList.reserve(Ybus.nonZeros());
-    for (int k=0; k < nb_bus_solver; ++k){
-        if(ybus_to_me(k) == -1) continue;  // I don't add anything to the slack bus
-        for (Eigen::SparseMatrix<cplx_type>::InnerIterator it(Ybus, k); it; ++it)
-        {
-            int row_res = static_cast<int>(it.row());  // TODO Eigen::Index here ?
-            if(ybus_to_me(row_res) == -1) continue;
-            row_res = ybus_to_me(row_res);
-            int col_res = static_cast<int>(it.col());  // should be k   // TODO Eigen::Index here ?
-            col_res = ybus_to_me(col_res);
-            tripletList.push_back(Eigen::Triplet<real_type> (row_res, col_res, std::real(it.value())));
-        }
-    }
-    dcYbus.setFromTriplets(tripletList.begin(), tripletList.end());
-    dcYbus.makeCompressed();
+    fill_dcYbus_noslack(nb_bus_solver, Ybus);
     
     #ifdef __COUT_TIMES
         std::cout << "\t dc: preproc: " << 1000. * timer_preproc.duration() << "ms" << std::endl;
@@ -85,7 +62,7 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     #endif // __COUT_TIMES
     bool just_factorize = false;
     if(need_factorize_){
-        ErrorType status_init = _linear_solver.initialize(dcYbus);
+        ErrorType status_init = _linear_solver.initialize(dcYbus_noslack_);
         if(status_init != ErrorType::NoError){
             err_ = status_init;
             return false;
@@ -95,17 +72,16 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     }
 
     // remove the slack bus from Sbus
-    RealVect dcSbus = RealVect::Constant(nb_bus_solver - slack_bus_ids_solver.size(), my_zero_);
+    dcSbus_noslack_ = RealVect::Constant(nb_bus_solver - slack_buses_ids_solver_.size(), my_zero_);
     for (int k=0; k < nb_bus_solver; ++k){
-        if(ybus_to_me(k) == -1) continue;  // I don't add anything to the slack bus
-        const int col_res = ybus_to_me(k);
-        dcSbus(col_res) = std::real(Sbus_tmp(k));
+        if(mat_bus_id_(k) == -1) continue;  // I don't add anything to the slack bus
+        const int col_res = mat_bus_id_(k);
+        dcSbus_noslack_(col_res) = std::real(Sbus_tmp(k));
     }
 
-    // solve for theta: Sbus = dcY . theta
-    RealVect Va_dc_without_slack = dcSbus;
-
-    ErrorType error = _linear_solver.solve(dcYbus, Va_dc_without_slack, just_factorize);
+    // solve for theta: Sbus = dcY . theta (make a copy to keep dcSbus_noslack_)
+    RealVect Va_dc_without_slack = dcSbus_noslack_;
+    ErrorType error = _linear_solver.solve(dcYbus_noslack_, Va_dc_without_slack, just_factorize);
     if(error != ErrorType::NoError){
         err_ = error;
         timer_total_nr_ += timer.duration();
@@ -135,18 +111,18 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     RealVect Va_dc = RealVect::Constant(nb_bus_solver, my_zero_);
     // fill Va from dc approx
     for (int ybus_id=0; ybus_id < nb_bus_solver; ++ybus_id){
-        if(ybus_to_me(ybus_id) == -1) continue;  // slack bus is handled elsewhere
-        const int bus_me = ybus_to_me(ybus_id);
+        if(mat_bus_id_(ybus_id) == -1) continue;  // slack bus is handled elsewhere
+        const int bus_me = mat_bus_id_(ybus_id);
         Va_dc(ybus_id) = Va_dc_without_slack(bus_me);
     }
-    Va_dc.array() += std::arg(V(slack_bus_ids_solver(0)));
+    Va_dc.array() += std::arg(V(slack_buses_ids_solver_(0)));
 
     // save the results
     Va_ = Va_dc;
 
     // add the Voltage setpoints of the generator
     Vm_ = V.array().abs();
-    Vm_(slack_bus_ids_solver) = V(slack_bus_ids_solver).array().abs();
+    Vm_(slack_buses_ids_solver_) = V(slack_buses_ids_solver_).array().abs();
 
     // now compute the resulting complex voltage
     V_ = (Va_.array().cos().template cast<cplx_type>() + my_i * Va_.array().sin().template cast<cplx_type>());
@@ -163,11 +139,57 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     return true;
 }
 
+template<class LinearSolver>
+void BaseDCSolver<LinearSolver>::fill_mat_bus_id(int nb_bus_solver){
+    mat_bus_id_ = Eigen::VectorXi::Constant(nb_bus_solver, -1);
+    // Eigen::VectorXi me_to_ybus = Eigen::VectorXi::Constant(nb_bus_solver - slack_bus_ids_solver.size(), -1);
+    int solver_id = 0;
+    for (int ybus_id=0; ybus_id < nb_bus_solver; ++ybus_id){
+        if(isin(ybus_id, slack_buses_ids_solver_)) continue;  // I don't add anything to the slack bus
+        mat_bus_id_(ybus_id) = solver_id;
+        // me_to_ybus(solver_id) = ybus_id;
+        ++solver_id;
+    }
+}
+
+template<class LinearSolver>
+void BaseDCSolver<LinearSolver>::fill_dcYbus_noslack(int nb_bus_solver, const Eigen::SparseMatrix<cplx_type> & ref_mat){
+    // TODO see if "prune" might work here https://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html#title29
+    remove_slack_buses(nb_bus_solver, ref_mat, dcYbus_noslack_);
+}
+
+template<class LinearSolver>
+template<typename ref_mat_type>  // ref_mat_type should be `real_type` or `cplx_type`
+void BaseDCSolver<LinearSolver>::remove_slack_buses(int nb_bus_solver, const Eigen::SparseMatrix<ref_mat_type> & ref_mat, Eigen::SparseMatrix<real_type> & res_mat){
+    res_mat = Eigen::SparseMatrix<real_type>(nb_bus_solver - 1, nb_bus_solver - 1);  // TODO dist slack: -1 or -mat_bus_id_.size() here ????
+    std::vector<Eigen::Triplet<real_type> > tripletList;
+    tripletList.reserve(ref_mat.nonZeros());
+    for (int k=0; k < nb_bus_solver; ++k){
+        if(mat_bus_id_(k) == -1) continue;  // I don't add anything to the slack bus
+        for (typename Eigen::SparseMatrix<ref_mat_type>::InnerIterator it(ref_mat, k); it; ++it)
+        {
+            int row_res = static_cast<int>(it.row());  // TODO Eigen::Index here ?
+            if(mat_bus_id_(row_res) == -1) continue;
+            row_res = mat_bus_id_(row_res);
+            int col_res = static_cast<int>(it.col());  // should be k   // TODO Eigen::Index here ?
+            col_res = mat_bus_id_(col_res);
+            tripletList.push_back(Eigen::Triplet<real_type> (row_res, col_res, std::real(it.value())));
+        }
+    }
+    res_mat.setFromTriplets(tripletList.begin(), tripletList.end());
+    res_mat.makeCompressed();
+}
+
 template<class LinearSolver>
 void BaseDCSolver<LinearSolver>::reset(){
     BaseSolver::reset();
     _linear_solver.reset();
     need_factorize_ = true;
+    dcSbus_noslack_ = RealVect();
+    dcYbus_noslack_ = Eigen::SparseMatrix<real_type>();
+    my_pv_ = Eigen::VectorXi();
+    slack_buses_ids_solver_ = Eigen::VectorXi();
+    mat_bus_id_ = Eigen::VectorXi();
 }
 
 template<class LinearSolver>
@@ -175,19 +197,19 @@ Eigen::SparseMatrix<real_type> BaseDCSolver<LinearSolver>::get_ptdf(){
     // TODO
     // Bf (nb_branch, nb_bus) : en dc un truc du genre 1 / x / tap for (1..nb_branch, from_bus)
     // and -1. / x / tap for (1..nb_branch, to_bus) 
-    return Ybus_noslack_;
+    return dcYbus_noslack_;
 }
 
 template<class LinearSolver>
 Eigen::SparseMatrix<real_type> BaseDCSolver<LinearSolver>::get_lodf(){
     // TODO
-    return Ybus_noslack_;
+    return dcYbus_noslack_;
 
 }
 
 template<class LinearSolver>
 Eigen::SparseMatrix<real_type> BaseDCSolver<LinearSolver>::get_bsdf(){
     // TODO
-    return Ybus_noslack_;
+    return dcYbus_noslack_;
 
 }
diff --git a/src/KLUSolver.cpp b/src/KLUSolver.cpp
index 316712a..fe8efc4 100644
--- a/src/KLUSolver.cpp
+++ b/src/KLUSolver.cpp
@@ -10,6 +10,8 @@
 
 #include <iostream>
 
+const bool KLULinearSolver::CAN_SOLVE_MAT = false;
+
 ErrorType KLULinearSolver::reset(){
     klu_free_symbolic(&symbolic_, &common_);
     klu_free_numeric(&numeric_, &common_);
diff --git a/src/KLUSolver.h b/src/KLUSolver.h
index c7aff02..2eefe06 100644
--- a/src/KLUSolver.h
+++ b/src/KLUSolver.h
@@ -48,6 +48,9 @@ class KLULinearSolver
         ErrorType initialize(Eigen::SparseMatrix<real_type>& J);
         ErrorType solve(Eigen::SparseMatrix<real_type>& J, RealVect & b, bool has_just_been_inialized);
 
+        // can this linear solver solve problem where RHS is a matrix
+        static const bool CAN_SOLVE_MAT;
+        
     private:
         // solver initialization
         klu_symbolic* symbolic_;
diff --git a/src/NICSLUSolver.cpp b/src/NICSLUSolver.cpp
index 18fed9e..357a672 100644
--- a/src/NICSLUSolver.cpp
+++ b/src/NICSLUSolver.cpp
@@ -9,6 +9,8 @@
 #include "NICSLUSolver.h"
 #include <iostream>
 
+const bool NICSLULinearSolver::CAN_SOLVE_MAT = false;
+
 ErrorType NICSLULinearSolver::reset(){
     // free everything
     solver_.Free();
diff --git a/src/NICSLUSolver.h b/src/NICSLUSolver.h
index 825bc98..481b4a0 100644
--- a/src/NICSLUSolver.h
+++ b/src/NICSLUSolver.h
@@ -56,6 +56,9 @@ class NICSLULinearSolver
         ErrorType initialize(Eigen::SparseMatrix<real_type> & J);
         ErrorType solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized);
 
+        // can this linear solver solve problem where RHS is a matrix
+        static const bool CAN_SOLVE_MAT;
+
         // prevent copy and assignment
         NICSLULinearSolver(const NICSLULinearSolver & other) = delete;
         NICSLULinearSolver & operator=( const NICSLULinearSolver & ) = delete;
diff --git a/src/SparseLUSolver.cpp b/src/SparseLUSolver.cpp
index 640a879..98c3149 100644
--- a/src/SparseLUSolver.cpp
+++ b/src/SparseLUSolver.cpp
@@ -10,6 +10,8 @@
 
 #include <iostream>
 
+const bool SparseLULinearSolver::CAN_SOLVE_MAT = true;
+
 ErrorType SparseLULinearSolver::initialize(const Eigen::SparseMatrix<real_type> & J){
     // default Eigen representation: column major, which is good for klu !
     // J is const here
diff --git a/src/SparseLUSolver.h b/src/SparseLUSolver.h
index 34849a5..768fd67 100644
--- a/src/SparseLUSolver.h
+++ b/src/SparseLUSolver.h
@@ -37,6 +37,8 @@ class SparseLULinearSolver
         ErrorType solve(const Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized);
         ErrorType reset(){ return ErrorType::NoError; }
 
+        // can this linear solver solve problem where RHS is a matrix
+        static const bool CAN_SOLVE_MAT;
     private:
         // solver initialization
         Eigen::SparseLU<Eigen::SparseMatrix<real_type>, Eigen::COLAMDOrdering<int> >  solver_;