diff --git a/DynamicModel.cc b/DynamicModel.cc
index a23253146e992e85ac6cf159d49be19a49fabca4..995a20286d7fefcf864e4024f4dbc51c34556a0a 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -2151,18 +2151,18 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
       computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
 
-      computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered, prologue, epilogue);
+      computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
 
       map<pair<int, pair<int, int> >, NodeID> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
 
-      equation_type_and_normalized_equation = equationTypeDetermination(equations, first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
+      equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
       cout << "Finding the optimal block decomposition of the model ...\n";
 
       if (prologue+epilogue < (unsigned int) equation_number())
-        computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, prologue, epilogue, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered);
+        computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered);
 
-      block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, prologue, epilogue, blocks, equations, equation_type_and_normalized_equation, variable_reordered, equation_reordered);
+      block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered);
 
       printBlockDecomposition(blocks);
 
diff --git a/ModelTree.cc b/ModelTree.cc
index 2e2c403e64a0093c1a0ca85c6c94187bef14f827..3fc5d76a432bba33d8a268d13fddfb6b879965a5 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -246,7 +246,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_type &eval_context, jaco
 {
   int nb_elements_contemparenous_Jacobian = 0;
   set<pair<int, int> > jacobian_elements_to_delete;
-  for (first_derivatives_type::iterator it = first_derivatives.begin();
+  for (first_derivatives_type::const_iterator it = first_derivatives.begin();
        it != first_derivatives.end(); it++)
     {
       int deriv_id = it->first.second;
@@ -304,7 +304,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_type &eval_context, jaco
 }
 
 void
-ModelTree::computePrologueAndEpilogue(jacob_map &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered, unsigned int &prologue, unsigned int &epilogue)
+ModelTree::computePrologueAndEpilogue(const jacob_map &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered)
 {
   vector<int> eq2endo(equation_number(), 0);
   equation_reordered.resize(equation_number());
@@ -411,7 +411,7 @@ ModelTree::computePrologueAndEpilogue(jacob_map &static_jacobian_arg, vector<int
 }
 
 t_equation_type_and_normalized_equation
-ModelTree::equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int mfs)
+ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
 {
   NodeID lhs, rhs;
   BinaryOpNode *eq_node;
@@ -419,14 +419,13 @@ ModelTree::equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair
   t_equation_type_and_normalized_equation V_Equation_Simulation_Type(equations.size());
   for (unsigned int i = 0; i < equations.size(); i++)
     {
-      temporary_terms.clear();
       int eq = Index_Equ_IM[i];
       int var = Index_Var_IM[i];
       eq_node = equations[eq];
       lhs = eq_node->get_arg1();
       rhs = eq_node->get_arg2();
       Equation_Simulation_Type = E_SOLVE;
-      map<pair<int, pair<int, int> >, NodeID>::iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
+      map<pair<int, pair<int, int> >, NodeID>::const_iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
       pair<bool, NodeID> res;
       if (derivative != first_order_endo_derivatives.end())
         {
@@ -460,7 +459,7 @@ ModelTree::equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair
 }
 
 void
-ModelTree::getVariableLeadLagByBlock(dynamic_jacob_map &dynamic_jacobian, vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_lag_lead_vector &equation_lead_lag, t_lag_lead_vector &variable_lead_lag, vector<int> equation_reordered, vector<int> variable_reordered) const
+ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, t_lag_lead_vector &equation_lead_lag, t_lag_lead_vector &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const
 {
   int nb_endo = symbol_table.endo_nbr();
   variable_lead_lag = t_lag_lead_vector(nb_endo, make_pair(0, 0));
@@ -504,7 +503,7 @@ ModelTree::getVariableLeadLagByBlock(dynamic_jacob_map &dynamic_jacobian, vector
 }
 
 void
-ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian, int prologue, int epilogue, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, t_equation_type_and_normalized_equation &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const
+ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map &static_jacobian, const dynamic_jacob_map &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, const t_equation_type_and_normalized_equation &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const
 {
   int nb_var = variable_reordered.size();
   int n = nb_var - prologue - epilogue;
@@ -573,7 +572,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(jacob_map &
 
   t_lag_lead_vector equation_lag_lead, variable_lag_lead;
 
-  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, prologue, epilogue, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
+  getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
 
   vector<int> tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
   int order = prologue;
@@ -637,7 +636,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(jacob_map &
 }
 
 void
-ModelTree::printBlockDecomposition(vector<pair<int, int> > blocks)
+ModelTree::printBlockDecomposition(const vector<pair<int, int> > &blocks) const
 {
   int largest_block = 0;
   int Nb_SimulBlocks = 0;
@@ -666,7 +665,7 @@ ModelTree::printBlockDecomposition(vector<pair<int, int> > blocks)
 }
 
 t_block_type_firstequation_size_mfs
-ModelTree::reduceBlocksAndTypeDetermination(dynamic_jacob_map &dynamic_jacobian, int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_equation_type_and_normalized_equation &Equation_Type, vector<int> &variable_reordered, vector<int> &equation_reordered)
+ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map &dynamic_jacobian, const vector<pair<int, int> > &blocks, const t_equation_type_and_normalized_equation &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered)
 {
   int i = 0;
   int count_equ = 0, blck_count_simult = 0;
@@ -786,7 +785,7 @@ ModelTree::reduceBlocksAndTypeDetermination(dynamic_jacob_map &dynamic_jacobian,
 }
 
 vector<bool>
-ModelTree::BlockLinear(t_blocks_derivatives &blocks_derivatives, vector<int> &variable_reordered)
+ModelTree::BlockLinear(const t_blocks_derivatives &blocks_derivatives, const vector<int> &variable_reordered) const
 {
   unsigned int nb_blocks = getNbBlocks();
   vector<bool> blocks_linear(nb_blocks, true);
diff --git a/ModelTree.hh b/ModelTree.hh
index 094f4666f5d115499f382a417f61bf6dc18aeb54..5415bf9cca65f1fdf490933a8cb5ef76bb729049 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -163,19 +163,19 @@ protected:
   //! Evaluate the jacobian and suppress all the elements below the cutoff
   void evaluateAndReduceJacobian(const eval_context_type &eval_context, jacob_map &contemporaneous_jacobian, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian, double cutoff, bool verbose);
   //! Search the equations and variables belonging to the prologue and the epilogue of the model
-  void computePrologueAndEpilogue(jacob_map &static_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, unsigned int &prologue, unsigned int &epilogue);
+  void computePrologueAndEpilogue(const jacob_map &static_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered);
   //! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
-  t_equation_type_and_normalized_equation equationTypeDetermination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int mfs);
+  t_equation_type_and_normalized_equation equationTypeDetermination(const map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
   //! Compute the block decomposition and for a non-recusive block find the minimum feedback set
-  void computeBlockDecompositionAndFeedbackVariablesForEachBlock(jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, t_equation_type_and_normalized_equation &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const;
+  void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map &static_jacobian, const dynamic_jacob_map &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, const t_equation_type_and_normalized_equation &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const;
   //! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
-  t_block_type_firstequation_size_mfs reduceBlocksAndTypeDetermination(dynamic_jacob_map &dynamic_jacobian, int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_equation_type_and_normalized_equation &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
+  t_block_type_firstequation_size_mfs reduceBlocksAndTypeDetermination(const dynamic_jacob_map &dynamic_jacobian, const vector<pair<int, int> > &blocks, const t_equation_type_and_normalized_equation &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered);
   //! Determine the maximum number of lead and lag for the endogenous variable in a bloc
-  void getVariableLeadLagByBlock(dynamic_jacob_map &dynamic_jacobian, vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue, t_lag_lead_vector &equation_lead_lag, t_lag_lead_vector &variable_lead_lag, vector<int> equation_reordered, vector<int> variable_reordered) const;
+  void getVariableLeadLagByBlock(const dynamic_jacob_map &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, t_lag_lead_vector &equation_lead_lag, t_lag_lead_vector &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const;
   //! Print an abstract of the block structure of the model
-  void printBlockDecomposition(vector<pair<int, int> > blocks);
+  void printBlockDecomposition(const vector<pair<int, int> > &blocks) const;
   //! Determine for each block if it is linear or not
-  vector<bool> BlockLinear(t_blocks_derivatives &blocks_derivatives, vector<int> &variable_reordered);
+  vector<bool> BlockLinear(const t_blocks_derivatives &blocks_derivatives, const vector<int> &variable_reordered) const;
 
   virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException) = 0;
   virtual int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException) = 0;
diff --git a/StaticModel.cc b/StaticModel.cc
index 5df774e88d9c456466b0ac47e0aa30964f133533..64935f7f26247767fe9686eff4a91d9b6df6debe 100644
--- a/StaticModel.cc
+++ b/StaticModel.cc
@@ -866,18 +866,18 @@ StaticModel::computingPass(const eval_context_type &eval_context, bool no_tmp_te
 
       computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
 
-      computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered, prologue, epilogue);
+      computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
 
       map<pair<int, pair<int, int> >, NodeID> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
 
-      equation_type_and_normalized_equation = equationTypeDetermination(equations, first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
+      equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
       cout << "Finding the optimal block decomposition of the model ...\n";
 
       if (prologue+epilogue < (unsigned int) equation_number())
-        computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, prologue, epilogue, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, false, mfs, inv_equation_reordered, inv_variable_reordered);
+        computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, false, mfs, inv_equation_reordered, inv_variable_reordered);
 
-      block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, prologue, epilogue, blocks, equations, equation_type_and_normalized_equation, variable_reordered, equation_reordered);
+      block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered);
 
       printBlockDecomposition(blocks);