diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index e0308aab..a4ae1b2e 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1249,7 +1249,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
                                block_size,
                                endo_idx_block2orig,
                                eq_idx_block2orig,
-                               blocks_linear[block],
+                               blocks[block].linear,
                                symbol_table.endo_nbr(),
                                block_max_lag,
                                block_max_lead,
@@ -2083,7 +2083,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                             << "  end;" << endl
                             << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
                             << ", y, x, params, steady_state, y_index, " << nze
-                            << ", options_.periods, " << blocks_linear[block]
+                            << ", options_.periods, " << blocks[block].linear
                             << ", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0, M_, options_, oo_);" << endl
                             << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
                             << "  if any(isnan(tmp) | isinf(tmp))" << endl
@@ -2115,7 +2115,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                             << "  end;" << endl
                             << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
                             <<", y, x, params, steady_state, y_index, " << nze
-                            <<", options_.periods, " << blocks_linear[block]
+                            <<", options_.periods, " << blocks[block].linear
                             <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0, M_, options_, oo_);" << endl
                             << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
                             << "  if any(isnan(tmp) | isinf(tmp))" << endl
@@ -2147,7 +2147,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                             <<", y, x, params, steady_state, y_index, " << nze
                             <<", options_.periods, " << max_leadlag_block[block].first
                             <<", " << max_leadlag_block[block].second
-                            <<", " << blocks_linear[block]
+                            <<", " << blocks[block].linear
                             <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, options_, M_, oo_);" << endl
                             << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
                             << "  if any(isnan(tmp) | isinf(tmp))" << endl
@@ -4793,30 +4793,22 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
       computeParamsDerivatives(paramsDerivsOrder);
     }
 
-  jacob_map_t contemporaneous_jacobian, static_jacobian;
-  map<tuple<int, int, int>, expr_t> first_order_endo_derivatives;
-  // For each simultaneous block contains pair<Size, Num_Feedback_variable>
-  vector<pair<int, int>> simblock_size;
-  vector<int> n_static, n_forward, n_backward, n_mixed;
 
   if (linear_decomposition)
     {
-      first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
-      is_equation_linear = equationLinear(first_order_endo_derivatives);
+      auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
+      equationLinear(first_order_endo_derivatives);
 
-      tie(contemporaneous_jacobian, static_jacobian) = evaluateAndReduceJacobian(eval_context, cutoff, false);
+      auto [contemporaneous_jacobian, static_jacobian] = evaluateAndReduceJacobian(eval_context, cutoff, false);
 
       if (!computeNaturalNormalization())
         computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
 
-      tie(simblock_size, n_static, n_forward, n_backward, n_mixed)
-        = select_non_linear_equations_and_variables(is_equation_linear);
+      select_non_linear_equations_and_variables();
 
       equationTypeDetermination(first_order_endo_derivatives, 0);
-      prologue = 0;
-      epilogue = 0;
 
-      reduceBlocksAndTypeDetermination(simblock_size, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
+      reduceBlocksAndTypeDetermination(linear_decomposition);
 
       computeChainRuleJacobian();
 
@@ -4830,26 +4822,23 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
       if (!no_tmp_terms)
         computeTemporaryTermsOrdered();
     }
-
-  if (block)
+  else if (block)
     {
-      tie(contemporaneous_jacobian, static_jacobian) = evaluateAndReduceJacobian(eval_context, cutoff, false);
+      auto [contemporaneous_jacobian, static_jacobian] = evaluateAndReduceJacobian(eval_context, cutoff, false);
 
       computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian);
 
       computePrologueAndEpilogue(static_jacobian);
 
-      first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
+      auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
 
       equationTypeDetermination(first_order_endo_derivatives, mfs);
 
       cout << "Finding the optimal block decomposition of the model ..." << endl;
 
-      lag_lead_vector_t variable_lag_lead;
+      auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, false);
 
-      tie(simblock_size, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false);
-
-      reduceBlocksAndTypeDetermination(simblock_size, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, linear_decomposition);
+      reduceBlocksAndTypeDetermination(linear_decomposition);
 
       printBlockDecomposition();
 
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index db2c35c4..d9e13f0d 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -162,19 +162,16 @@ ModelTree::ModelTree(const ModelTree &m) :
   eq_idx_orig2block{m.eq_idx_orig2block},
   endo_idx_orig2block{m.endo_idx_orig2block},
   map_idx{m.map_idx},
-  block_type_firstequation_size_mfs{m.block_type_firstequation_size_mfs},
-  blocks_linear{m.blocks_linear},
-  block_col_type{m.block_col_type},
   endo_max_leadlag_block{m.endo_max_leadlag_block},
   other_endo_max_leadlag_block{m.other_endo_max_leadlag_block},
   exo_max_leadlag_block{m.exo_max_leadlag_block},
   exo_det_max_leadlag_block{m.exo_det_max_leadlag_block},
   max_leadlag_block{m.max_leadlag_block},
+  blocks{m.blocks},
   is_equation_linear{m.is_equation_linear},
   endo2eq{m.endo2eq},
   epilogue{m.epilogue},
   prologue{m.prologue},
-  block_lag_lead{m.block_lag_lead},
   cutoff{m.cutoff},
   mfs{m.mfs}
 {
@@ -215,25 +212,21 @@ ModelTree::operator=(const ModelTree &m)
   first_chain_rule_derivatives.clear();
   map_idx = m.map_idx;
   equation_type_and_normalized_equation.clear();
-  block_type_firstequation_size_mfs = m.block_type_firstequation_size_mfs;
   blocks_derivatives.clear();
-  blocks_linear = m.blocks_linear;
   derivative_endo.clear();
   derivative_other_endo.clear();
   derivative_exo.clear();
   derivative_exo_det.clear();
-  block_col_type = m.block_col_type;
   endo_max_leadlag_block = m.endo_max_leadlag_block;
   other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
   exo_max_leadlag_block = m.exo_max_leadlag_block;
   exo_det_max_leadlag_block = m.exo_det_max_leadlag_block;
   max_leadlag_block = m.max_leadlag_block;
-
+  blocks = m.blocks;
   is_equation_linear = m.is_equation_linear;
   endo2eq = m.endo2eq;
   epilogue = m.epilogue;
   prologue = m.prologue;
-  block_lag_lead = m.block_lag_lead;
   cutoff = m.cutoff;
   mfs = m.mfs;
 
@@ -456,43 +449,46 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, double
   return { contemporaneous_jacobian, static_jacobian };
 }
 
-tuple<vector<pair<int, int>>, vector<int>, vector<int>, vector<int>, vector<int>>
-ModelTree::select_non_linear_equations_and_variables(const vector<bool> &is_equation_linear)
+void
+ModelTree::select_non_linear_equations_and_variables()
 {
-  vector<int> eq2endo(equations.size(), 0);
-  int num = 0;
-  for (auto it : endo2eq)
-    if (!is_equation_linear[it])
-      num++;
-  vector<int> endo2block(endo2eq.size(), 1);
-  int i = 0, j = 0;
-  for (auto it : endo2eq)
-    if (!is_equation_linear[it])
-      {
-        eq_idx_block2orig[i] = it;
-        endo_idx_block2orig[i] = j;
-        endo2block[j] = 0;
-        i++;
-        j++;
-      }
-  auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block);
-  vector<int> n_static(endo2eq.size(), 0), n_forward(endo2eq.size(), 0),
-    n_backward(endo2eq.size(), 0), n_mixed(endo2eq.size(), 0);
-  for (int i = 0; i < static_cast<int>(endo2eq.size()); i++)
+  prologue = 0;
+  epilogue = 0;
+
+  vector<int> endo2block(endo2eq.size(), 1); // The 1 is a dummy value, distinct from 0
+  int i = 0;
+  for (int endo = 0; endo < static_cast<int>(endo2eq.size()); endo++)
     {
-      if (variable_lag_lead[endo_idx_block2orig[i]].first != 0 && variable_lag_lead[endo_idx_block2orig[i]].second != 0)
-        n_mixed[i]++;
-      else if (variable_lag_lead[endo_idx_block2orig[i]].first == 0 && variable_lag_lead[endo_idx_block2orig[i]].second != 0)
-        n_forward[i]++;
-      else if (variable_lag_lead[endo_idx_block2orig[i]].first != 0 && variable_lag_lead[endo_idx_block2orig[i]].second == 0)
-        n_backward[i]++;
-      else if (variable_lag_lead[endo_idx_block2orig[i]].first == 0 && variable_lag_lead[endo_idx_block2orig[i]].second == 0)
-        n_static[i]++;
+      int eq = endo2eq[endo];
+      if (!is_equation_linear[eq])
+        {
+          eq_idx_block2orig[i] = eq;
+          endo_idx_block2orig[i] = endo;
+          endo2block[i] = 0;
+          i++;
+        }
     }
-  cout.flush();
-  vector<pair<int, int>> simblock_size(1, {i, i});
   updateReverseVariableEquationOrderings();
-  return { simblock_size, n_static, n_forward, n_backward, n_mixed };
+
+  blocks.clear();
+  blocks.resize(1);
+  blocks[0].size = i;
+  blocks[0].mfs_size = i;
+
+  auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block);
+
+  for (int i = 0; i < blocks[0].size; i++)
+    {
+      auto [max_lag, max_lead] = variable_lag_lead[endo_idx_block2orig[i]];
+      if (max_lag != 0 && max_lead != 0)
+        blocks[0].n_mixed++;
+      else if (max_lag == 0 && max_lead != 0)
+        blocks[0].n_forward++;
+      else if (max_lag != 0 && max_lead == 0)
+        blocks[0].n_backward++;
+      else
+        blocks[0].n_static++;
+    }
 }
 
 bool
@@ -701,9 +697,8 @@ ModelTree::getVariableLeadLagByBlock(const vector<int> &endo2simblock) const
   return { equation_lag_lead, variable_lag_lead };
 }
 
-tuple<vector<pair<int, int>>, lag_lead_vector_t,
-      vector<int>, vector<int>, vector<int>, vector<int>>
-ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_)
+lag_lead_vector_t
+ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, bool verbose_)
 {
   int nb_var = symbol_table.endo_nbr();
   int nb_simvars = nb_var - prologue - epilogue;
@@ -731,13 +726,28 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
      index, starting from 0, in recursive order */
   auto [num_simblocks, endo2simblock] = G.sortedStronglyConnectedComponents();
 
-  vector<pair<int, int>> simblock_size(num_simblocks, { 0, 0 });
+  int num_blocks = prologue+num_simblocks+epilogue;
 
-  // equations belonging to the block
+  blocks.clear();
+  blocks.resize(num_blocks);
+
+  // Initialize size and mfs_size for prologue and epilogue
+  for (int i = 0; i < prologue; i++)
+    {
+      blocks[i].size = 1;
+      blocks[i].mfs_size = 1;
+    }
+  for (int i = 0; i < epilogue; i++)
+    {
+      blocks[prologue+num_simblocks+i].size = 1;
+      blocks[prologue+num_simblocks+i].mfs_size = 1;
+    }
+
+  // Compute size and list of equations for simultaneous blocks
   vector<set<int>> eqs_in_simblock(num_simblocks);
   for (int i = 0; i < static_cast<int>(endo2simblock.size()); i++)
     {
-      simblock_size[endo2simblock[i]].first++;
+      blocks[prologue+endo2simblock[i]].size++;
       eqs_in_simblock[endo2simblock[i]].insert(i);
     }
 
@@ -746,7 +756,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
   /* Add a loop on vertices which could not be normalized or vertices related
      to lead/lag variables. This forces those vertices to belong to the feedback set */
   for (int i = 0; i < nb_simvars; i++)
-    if (Equation_Type[eq_idx_block2orig[i+prologue]].first == EquationType::solve
+    if (equation_type_and_normalized_equation[eq_idx_block2orig[i+prologue]].first == EquationType::solve
         || variable_lag_lead[endo_idx_block2orig[i+prologue]].first > 0
         || variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0
         || equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0
@@ -754,23 +764,19 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
         || mfs == 0)
       add_edge(vertex(i, G), vertex(i, G), G);
 
-  int num_blocks = prologue+num_simblocks+epilogue;
   // Determines the dynamic structure of each equation
-  vector<int> n_static(num_blocks, 0), n_forward(num_blocks, 0),
-    n_backward(num_blocks, 0), n_mixed(num_blocks, 0);
-
   const vector<int> old_eq_idx_block2orig(eq_idx_block2orig), old_endo_idx_block2orig(endo_idx_block2orig);
   for (int i = 0; i < prologue; i++)
     {
       auto [max_lag, max_lead] = variable_lag_lead[old_endo_idx_block2orig[i]];
       if (max_lag != 0 && max_lead != 0)
-        n_mixed[i]++;
+        blocks[i].n_mixed++;
       else if (max_lag == 0 && max_lead != 0)
-        n_forward[i]++;
+        blocks[i].n_forward++;
       else if (max_lag != 0 && max_lead == 0)
-        n_backward[i]++;
+        blocks[i].n_backward++;
       else
-        n_static[i]++;
+        blocks[i].n_static++;
     }
 
   /* For each block, the minimum set of feedback variable is computed and the
@@ -783,7 +789,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
       auto subG = G.extractSubgraph(eqs_in_simblock[i]);
       auto feed_back_vertices = subG.minimalSetOfFeedbackVertices();
       auto v_index1 = get(boost::vertex_index1, subG);
-      simblock_size[i].second = feed_back_vertices.size();
+      blocks[prologue+i].mfs_size = feed_back_vertices.size();
       auto reordered_vertices = subG.reorderRecursiveVariables(feed_back_vertices);
 
       // First we have the recursive equations conditional on feedback variables
@@ -799,22 +805,22 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
                            };
             if (j == 2 && max_lag != 0 && max_lead != 0)
               {
-                n_mixed[prologue+i]++;
+                blocks[prologue+i].n_mixed++;
                 reorder();
               }
             else if (j == 3 && max_lag == 0 && max_lead != 0)
               {
-                n_forward[prologue+i]++;
+                blocks[prologue+i].n_forward++;
                 reorder();
               }
             else if (j == 1 && max_lag != 0 && max_lead == 0)
               {
-                n_backward[prologue+i]++;
+                blocks[prologue+i].n_backward++;
                 reorder();
               }
             else if (j == 0 && max_lag == 0 && max_lead == 0)
               {
-                n_static[prologue+i]++;
+                blocks[prologue+i].n_static++;
                 reorder();
               }
           }
@@ -833,22 +839,22 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
                            };
             if (j == 2 && max_lag != 0 && max_lead != 0)
               {
-                n_mixed[prologue+i]++;
+                blocks[prologue+i].n_mixed++;
                 reorder();
               }
             else if (j == 3 && max_lag == 0 && max_lead != 0)
               {
-                n_forward[prologue+i]++;
+                blocks[prologue+i].n_forward++;
                 reorder();
               }
             else if (j == 1 && max_lag != 0 && max_lead == 0)
               {
-                n_backward[prologue+i]++;
+                blocks[prologue+i].n_backward++;
                 reorder();
               }
             else if (j == 0 && max_lag == 0 && max_lead == 0)
               {
-                n_static[prologue+i]++;
+                blocks[prologue+i].n_static++;
                 reorder();
               }
           }
@@ -858,18 +864,18 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
     {
       auto [max_lag, max_lead] = variable_lag_lead[old_endo_idx_block2orig[prologue+nb_simvars+i]];
       if (max_lag != 0 && max_lead != 0)
-        n_mixed[prologue+num_simblocks+i]++;
+        blocks[prologue+num_simblocks+i].n_mixed++;
       else if (max_lag == 0 && max_lead != 0)
-        n_forward[prologue+num_simblocks+i]++;
+        blocks[prologue+num_simblocks+i].n_forward++;
       else if (max_lag != 0 && max_lead == 0)
-        n_backward[prologue+num_simblocks+i]++;
+        blocks[prologue+num_simblocks+i].n_backward++;
       else
-        n_static[prologue+num_simblocks+i]++;
+        blocks[prologue+num_simblocks+i].n_static++;
     }
 
   updateReverseVariableEquationOrderings();
 
-  return { simblock_size, variable_lag_lead, n_static, n_forward, n_backward, n_mixed };
+  return variable_lag_lead;
 }
 
 void
@@ -904,185 +910,152 @@ ModelTree::printBlockDecomposition() const
 }
 
 void
-ModelTree::reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &simblock_size, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &n_static, const vector<int> &n_forward, const vector<int> &n_backward, const vector<int> &n_mixed, bool linear_decomposition)
+ModelTree::reduceBlocksAndTypeDetermination(bool linear_decomposition)
 {
-  int count_equ = 0, simblock_counter = 0;
-  block_type_firstequation_size_mfs.clear();
-  BlockSimulationType prev_Type = BlockSimulationType::unknown;
-  int eq = 0;
-  int num_simblocks = simblock_size.size();
-  for (int i = 0; i < prologue+num_simblocks+epilogue; i++)
+  for (int i = 0; i < static_cast<int>(blocks.size()); i++)
     {
-      int Blck_Size, MFS_Size;
-      int first_count_equ = count_equ;
-      if (i < prologue)
-        {
-          Blck_Size = 1;
-          MFS_Size = 1;
-        }
-      else if (i < prologue+num_simblocks)
-        {
-          Blck_Size = simblock_size[simblock_counter].first;
-          MFS_Size = simblock_size[simblock_counter].second;
-          simblock_counter++;
-        }
-      else if (i < prologue+num_simblocks+epilogue)
-        {
-          Blck_Size = 1;
-          MFS_Size = 1;
-        }
+      int first_eq = (i == 0) ? 0 : blocks[i-1].first_equation+blocks[i-1].size;
 
-      int Lag = 0, Lead = 0;
-      for (count_equ = first_count_equ; count_equ < Blck_Size+first_count_equ; count_equ++)
+      /* Compute the maximum lead and lag across all endogenous that appear in
+         this block and that belong to it */
+      int max_lag = 0, max_lead = 0;
+      for (int eq = first_eq; eq < first_eq+blocks[i].size; eq++)
         {
           set<pair<int, int>> endos_and_lags;
-          equations[eq_idx_block2orig[count_equ]]->collectEndogenous(endos_and_lags);
-          for (const auto &[curr_variable, curr_lag] : endos_and_lags)
+          equations[eq_idx_block2orig[eq]]->collectEndogenous(endos_and_lags);
+          for (const auto &[endo, lag] : endos_and_lags)
             {
               if (linear_decomposition)
                 {
-                  if (dynamic_jacobian.find({ curr_lag, eq_idx_block2orig[count_equ], curr_variable }) != dynamic_jacobian.end())
+                  if (dynamic_jacobian.find({ lag, eq_idx_block2orig[eq], endo })
+                      != dynamic_jacobian.end())
                     {
-                      if (curr_lag > Lead)
-                        Lead = curr_lag;
-                      else if (-curr_lag > Lag)
-                        Lag = -curr_lag;
+                      max_lead = max(lag, max_lead);
+                      max_lag = max(-lag, max_lag);
                     }
                 }
               else
                 {
-                  if (find(endo_idx_block2orig.begin()+first_count_equ, endo_idx_block2orig.begin()+(first_count_equ+Blck_Size), curr_variable)
-                      != endo_idx_block2orig.begin()+(first_count_equ+Blck_Size)
-                      && dynamic_jacobian.find({ curr_lag, eq_idx_block2orig[count_equ], curr_variable }) != dynamic_jacobian.end())
+                  if (find(endo_idx_block2orig.begin()+first_eq,
+                           endo_idx_block2orig.begin()+first_eq+blocks[i].size, endo)
+                      != endo_idx_block2orig.begin()+first_eq+blocks[i].size
+                      && dynamic_jacobian.find({ lag, eq_idx_block2orig[eq], endo })
+                      != dynamic_jacobian.end())
                     {
-                      if (curr_lag > Lead)
-                        Lead = curr_lag;
-                      else if (-curr_lag > Lag)
-                        Lag = -curr_lag;
+                      max_lead = max(lag, max_lead);
+                      max_lag = max(-lag, max_lag);
                     }
                 }
             }
         }
+
+      // Determine the block type
       BlockSimulationType Simulation_Type;
-      if (Lag > 0 && Lead > 0)
+      if (max_lag > 0 && max_lead > 0)
         {
-          if (Blck_Size == 1)
+          if (blocks[i].size == 1)
             Simulation_Type = BlockSimulationType::solveTwoBoundariesSimple;
           else
             Simulation_Type = BlockSimulationType::solveTwoBoundariesComplete;
         }
-      else if (Blck_Size > 1)
+      else if (blocks[i].size > 1)
         {
-          if (Lead > 0)
+          if (max_lead > 0)
             Simulation_Type = BlockSimulationType::solveBackwardComplete;
           else
             Simulation_Type = BlockSimulationType::solveForwardComplete;
         }
       else
         {
-          if (Lead > 0)
+          if (max_lead > 0)
             Simulation_Type = BlockSimulationType::solveBackwardSimple;
           else
             Simulation_Type = BlockSimulationType::solveForwardSimple;
         }
-      int l_n_static = n_static[i];
-      int l_n_forward = n_forward[i];
-      int l_n_backward = n_backward[i];
-      int l_n_mixed = n_mixed[i];
-      if (Blck_Size == 1)
+
+      if (blocks[i].size == 1)
         {
-          if (Equation_Type[eq_idx_block2orig[eq]].first == EquationType::evaluate
-              || Equation_Type[eq_idx_block2orig[eq]].first == EquationType::evaluate_s)
+          // Determine if the block can simply be evaluated
+          if (equation_type_and_normalized_equation[eq_idx_block2orig[first_eq]].first == EquationType::evaluate
+              || equation_type_and_normalized_equation[eq_idx_block2orig[first_eq]].first == EquationType::evaluate_s)
             {
               if (Simulation_Type == BlockSimulationType::solveBackwardSimple)
                 Simulation_Type = BlockSimulationType::evaluateBackward;
               else if (Simulation_Type == BlockSimulationType::solveForwardSimple)
                 Simulation_Type = BlockSimulationType::evaluateForward;
             }
+
+          /* Try to merge this block with the previous one.
+             This is only possible if the two blocks can simply be evaluated
+             (in the same direction), and if the merge does not break the
+             restrictions on leads/lags. */
           if (i > 0)
             {
               bool is_lead = false, is_lag = false;
-              int c_Size = get<2>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
-              int first_equation = get<1>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
-              if (c_Size > 0
-                  && ((prev_Type == BlockSimulationType::evaluateForward && Simulation_Type == BlockSimulationType::evaluateForward && !is_lead)
-                      || (prev_Type == BlockSimulationType::evaluateBackward && Simulation_Type == BlockSimulationType::evaluateBackward && !is_lag)))
+              for (int j = blocks[i-1].first_equation;
+                   j < blocks[i-1].first_equation+blocks[i-1].size; j++)
                 {
-                  for (int j = first_equation; j < first_equation+c_Size; j++)
-                    {
-                      if (dynamic_jacobian.find({ -1, eq_idx_block2orig[eq], endo_idx_block2orig[j] })
-                          != dynamic_jacobian.end())
-                        is_lag = true;
-                      if (dynamic_jacobian.find({ +1, eq_idx_block2orig[eq], endo_idx_block2orig[j] })
-                          != dynamic_jacobian.end())
-                        is_lead = true;
-                    }
+                  if (dynamic_jacobian.find({ -1, eq_idx_block2orig[first_eq], endo_idx_block2orig[j] })
+                      != dynamic_jacobian.end())
+                    is_lag = true;
+                  if (dynamic_jacobian.find({ +1, eq_idx_block2orig[first_eq], endo_idx_block2orig[j] })
+                      != dynamic_jacobian.end())
+                    is_lead = true;
                 }
-              if ((prev_Type == BlockSimulationType::evaluateForward && Simulation_Type == BlockSimulationType::evaluateForward && !is_lead)
-                  || (prev_Type == BlockSimulationType::evaluateBackward && Simulation_Type == BlockSimulationType::evaluateBackward && !is_lag))
+
+              BlockSimulationType prev_Type = blocks[i-1].simulation_type;
+              if ((prev_Type == BlockSimulationType::evaluateForward
+                   && Simulation_Type == BlockSimulationType::evaluateForward
+                   && !is_lead)
+                  || (prev_Type == BlockSimulationType::evaluateBackward
+                      && Simulation_Type == BlockSimulationType::evaluateBackward
+                      && !is_lag))
                 {
-                  //merge the current block with the previous one
-                  BlockSimulationType c_Type = get<0>(block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1]);
-                  c_Size++;
-                  block_type_firstequation_size_mfs[block_type_firstequation_size_mfs.size()-1] = { c_Type, first_equation, c_Size, c_Size };
-                  if (block_lag_lead[block_type_firstequation_size_mfs.size()-1].first > Lag)
-                    Lag = block_lag_lead[block_type_firstequation_size_mfs.size()-1].first;
-                  if (block_lag_lead[block_type_firstequation_size_mfs.size()-1].second > Lead)
-                    Lead = block_lag_lead[block_type_firstequation_size_mfs.size()-1].second;
-                  block_lag_lead[block_type_firstequation_size_mfs.size()-1] = { Lag, Lead };
-                  auto tmp = block_col_type[block_col_type.size()-1];
-                  block_col_type[block_col_type.size()-1] = { get<0>(tmp)+l_n_static, get<1>(tmp)+l_n_forward, get<2>(tmp)+l_n_backward, get<3>(tmp)+l_n_mixed };
-                }
-              else
-                {
-                  block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
-                  block_lag_lead.emplace_back(Lag, Lead);
-                  block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
+                  // Merge the current block into the previous one
+                  blocks[i-1].size++;
+                  blocks[i-1].mfs_size = blocks[i-1].size;
+                  blocks[i-1].max_lag = max(blocks[i-1].max_lag, max_lag);
+                  blocks[i-1].max_lead = max(blocks[i-1].max_lead, max_lead);
+                  blocks[i-1].n_static += blocks[i].n_static;
+                  blocks[i-1].n_forward += blocks[i].n_forward;
+                  blocks[i-1].n_backward += blocks[i].n_backward;
+                  blocks[i-1].n_mixed += blocks[i].n_mixed;
+                  blocks.erase(blocks.begin()+i);
+                  i--;
+                  continue;
                 }
             }
-          else
-            {
-              block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
-              block_lag_lead.emplace_back(Lag, Lead);
-              block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
-            }
         }
-      else
-        {
-          block_type_firstequation_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
-          block_lag_lead.emplace_back(Lag, Lead);
-          block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
-        }
-      prev_Type = Simulation_Type;
-      eq += Blck_Size;
+
+      blocks[i].simulation_type = Simulation_Type;
+      blocks[i].first_equation = first_eq;
+      blocks[i].max_lag = max_lag;
+      blocks[i].max_lead = max_lead;
     }
 }
 
-vector<bool>
-ModelTree::equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives) const
+void
+ModelTree::equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives)
 {
-  vector<bool> is_linear(symbol_table.endo_nbr(), true);
-  for (const auto &it : first_order_endo_derivatives)
+  is_equation_linear.clear();
+  is_equation_linear.resize(symbol_table.endo_nbr(), true);
+  for (const auto &[indices, expr] : first_order_endo_derivatives)
     {
-      expr_t Id = it.second;
       set<pair<int, int>> endogenous;
-      Id->collectEndogenous(endogenous);
+      expr->collectEndogenous(endogenous);
       if (endogenous.size() > 0)
         {
-          int eq = get<0>(it.first);
-          is_linear[eq] = false;
+          int eq = get<0>(indices);
+          is_equation_linear[eq] = false;
         }
     }
-  return is_linear;
 }
 
 void
 ModelTree::determineLinearBlocks()
 {
-  int nb_blocks = getNbBlocks();
-  blocks_linear.clear();
-  blocks_linear.resize(nb_blocks, true);
-  for (int block = 0; block < nb_blocks; block++)
+  // Note that field “linear” in class BlockInfo defaults to true
+  for (int block = 0; block < getNbBlocks(); block++)
     {
       BlockSimulationType simulation_type = getBlockSimulationType(block);
       int block_size = getBlockSize(block);
@@ -1100,7 +1073,7 @@ ModelTree::determineLinearBlocks()
                   for (int l = 0; l < block_size; l++)
                     if (endogenous.find({ endo_idx_block2orig[first_variable_position+l], 0 }) != endogenous.end())
                       {
-                        blocks_linear[block] = false;
+                        blocks[block].linear = false;
                         goto the_end;
                       }
               }
@@ -1115,7 +1088,7 @@ ModelTree::determineLinearBlocks()
               for (int l = 0; l < block_size; l++)
                 if (endogenous.find({ endo_idx_block2orig[first_variable_position+l], lag }) != endogenous.end())
                   {
-                    blocks_linear[block] = false;
+                    blocks[block].linear = false;
                     goto the_end;
                   }
           }
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index ea722a4c..17ad90a7 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -55,9 +55,6 @@ using equation_type_and_normalized_equation_t = vector<pair<EquationType, expr_t
 //! Vector describing variables: max_lag in the block, max_lead in the block
 using lag_lead_vector_t = vector<pair<int, int>>;
 
-//! for each block contains tuple<Simulation_Type, first_equation, Block_Size, Recursive_part_Size>
-using block_type_firstequation_size_mfs_t = vector<tuple<BlockSimulationType, int, int, int>>;
-
 //! for a block contains derivatives tuple<block_equation_number, block_variable_number, lead_lag, expr_t>
 using block_derivatives_equation_variable_laglead_nodeid_t = vector<tuple<int, int, int, expr_t>>;
 
@@ -190,27 +187,32 @@ protected:
   //! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
   equation_type_and_normalized_equation_t equation_type_and_normalized_equation;
 
-  //! For each block contains tuple<Simulation_Type, first_equation, Block_Size, Recursive_part_Size>
-  block_type_firstequation_size_mfs_t block_type_firstequation_size_mfs;
-
   //! for all blocks derivatives description
   blocks_derivatives_t blocks_derivatives;
 
-  //! Vector indicating if the block is linear in endogenous variable (true) or not (false)
-  vector<bool> blocks_linear;
-
   //! Map the derivatives for a block tuple<lag, eq, var>
   using derivative_t = map<tuple<int, int, int>, expr_t>;
   //! Vector of derivative for each blocks
   vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
 
-  //! for each block described the number of static, forward, backward and mixed variables in the block
-  /*! tuple<static, forward, backward, mixed> */
-  vector<tuple<int, int, int, int>> block_col_type;
-
   //!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
   vector<pair<int, int>> endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
 
+  class BlockInfo
+  {
+  public:
+    BlockSimulationType simulation_type;
+    int first_equation; // Stores a block-ordered equation ID
+    int size{0};
+    int mfs_size{0}; // Size of the minimal feedback set
+    bool linear{true}; // Whether the block is linear in endogenous variable
+    int n_static{0}, n_forward{0}, n_backward{0}, n_mixed{0};
+    int max_lag{0}, max_lead{0};
+  };
+
+  // Stores various informations on the blocks
+  vector<BlockInfo> blocks;
+
   //! the file containing the model and the derivatives code
   ofstream code_file;
 
@@ -270,9 +272,6 @@ protected:
   //! number of equation in the prologue and in the epilogue
   int epilogue, prologue;
 
-  //! for each block contains pair< max_lag, max_lead>
-  lag_lead_vector_t block_lag_lead;
-
   /* Compute a pseudo-Jacobian whose all elements are either zero or one,
      depending on whether the variable symbolically appears in the equation */
   jacob_map_t computeSymbolicJacobian() const;
@@ -302,17 +301,22 @@ protected:
     dynamic_jacobian. Elements below the cutoff are discarded. External functions are evaluated to 1. */
   pair<jacob_map_t, jacob_map_t> evaluateAndReduceJacobian(const eval_context_t &eval_context, double cutoff, bool verbose);
   //! Select and reorder the non linear equations of the model
-  /*! Returns a tuple (blocks, n_static, n_forward, n_backward, n_mixed) */
-  tuple<vector<pair<int, int>>, vector<int>, vector<int>, vector<int>, vector<int>> select_non_linear_equations_and_variables(const vector<bool> &is_equation_linear);
+  void select_non_linear_equations_and_variables();
   //! Search the equations and variables belonging to the prologue and the epilogue of the model
   void computePrologueAndEpilogue(const jacob_map_t &static_jacobian);
   //! Determine the type of each equation of model and try to normalize the unnormalized equation
   void equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs);
-  //! Compute the block decomposition and for a non-recusive block find the minimum feedback set
-  /*! Returns a tuple (blocks, variable_lag_lead, n_static, n_forward, n_backward, n_mixed) */
-  tuple<vector<pair<int, int>>, lag_lead_vector_t, vector<int>, vector<int>, vector<int>, vector<int>> computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_);
-  //! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
-  void reduceBlocksAndTypeDetermination(const vector<pair<int, int>> &simblock_size, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &n_static, const vector<int> &n_forward, const vector<int> &n_backward, const vector<int> &n_mixed, bool linear_decomposition);
+  /* Compute the block decomposition and for a non-recusive block find the minimum feedback set
+
+     Initializes the “blocks” structure, and fills the following fields: size,
+     mfs_size, n_static, n_forward, n_backward, n_mixed. */
+  lag_lead_vector_t computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, bool verbose_);
+  /* Reduce the number of block by merging the same type of equations in the
+     prologue and the epilogue, and determine the type of each block.
+
+     Fills the following fields of the “blocks” structure: simulation_type,
+     first_equation, max_lead, max_lag. */
+  void reduceBlocksAndTypeDetermination(bool linear_decomposition);
   /* The 1st output gives, for each equation (in original order) the (max_lag,
      max_lead) across all endogenous that appear in the equation and that
      belong to the same block (i.e. those endogenous are solved in the same
@@ -323,7 +327,7 @@ protected:
      which it belongs. */
   pair<lag_lead_vector_t, lag_lead_vector_t> getVariableLeadLagByBlock(const vector<int> &endo2simblock) const;
   //! For each equation determine if it is linear or not
-  vector<bool> equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives) const;
+  void equationLinear(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives);
   //! Print an abstract of the block structure of the model
   void printBlockDecomposition() const;
   //! Determine for each block if it is linear or not
@@ -337,25 +341,25 @@ protected:
   int
   getNbBlocks() const
   {
-    return block_type_firstequation_size_mfs.size();
+    return blocks.size();
   };
   //! Determine the simulation type of each block
   BlockSimulationType
   getBlockSimulationType(int block_number) const
   {
-    return get<0>(block_type_firstequation_size_mfs[block_number]);
+    return blocks[block_number].simulation_type;
   };
   //! Return the first equation number of a block
   int
   getBlockFirstEquation(int block_number) const
   {
-    return get<1>(block_type_firstequation_size_mfs[block_number]);
+    return blocks[block_number].first_equation;
   };
   //! Return the size of the block block_number
   int
   getBlockSize(int block_number) const
   {
-    return get<2>(block_type_firstequation_size_mfs[block_number]);
+    return blocks[block_number].size;
   };
   //! Return the number of exogenous variable in the block block_number
   virtual int getBlockExoSize(int block_number) const = 0;
@@ -365,61 +369,62 @@ protected:
   int
   getBlockMfs(int block_number) const
   {
-    return get<3>(block_type_firstequation_size_mfs[block_number]);
+    return blocks[block_number].mfs_size;
   };
   //! Return the maximum lag in a block
   int
   getBlockMaxLag(int block_number) const
   {
-    return block_lag_lead[block_number].first;
+    return blocks[block_number].max_lag;
   };
   //! Return the maximum lead in a block
   int
   getBlockMaxLead(int block_number) const
   {
-    return block_lag_lead[block_number].second;
+    return blocks[block_number].max_lead;
   };
   inline void
   setBlockLeadLag(int block, int max_lag, int max_lead)
   {
-    block_lag_lead[block] = { max_lag, max_lead };
+    blocks[block].max_lag = max_lag;
+    blocks[block].max_lead = max_lead;
   };
 
   //! Return the type of equation (equation_number) belonging to the block block_number
   EquationType
   getBlockEquationType(int block_number, int equation_number) const
   {
-    return equation_type_and_normalized_equation[eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first;
+    return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].first;
   };
   //! Return true if the equation has been normalized
   bool
   isBlockEquationRenormalized(int block_number, int equation_number) const
   {
-    return equation_type_and_normalized_equation[eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first == EquationType::evaluate_s;
+    return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].first == EquationType::evaluate_s;
   };
   //! Return the expr_t of the equation equation_number belonging to the block block_number
   expr_t
   getBlockEquationExpr(int block_number, int equation_number) const
   {
-    return equations[eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]];
+    return equations[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]];
   };
   //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
   expr_t
   getBlockEquationRenormalizedExpr(int block_number, int equation_number) const
   {
-    return equation_type_and_normalized_equation[eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].second;
+    return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].second;
   };
   //! Return the original number of equation equation_number belonging to the block block_number
   int
   getBlockEquationID(int block_number, int equation_number) const
   {
-    return eq_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number];
+    return eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number];
   };
   //! Return the original number of variable variable_number belonging to the block block_number
   int
   getBlockVariableID(int block_number, int variable_number) const
   {
-    return endo_idx_block2orig[get<1>(block_type_firstequation_size_mfs[block_number])+variable_number];
+    return endo_idx_block2orig[getBlockFirstEquation(block_number)+variable_number];
   };
   //! Return the original number of the exogenous variable varexo_number belonging to the block block_number
   virtual int getBlockVariableExoID(int block_number, int variable_number) const = 0;
@@ -427,13 +432,13 @@ protected:
   int
   getBlockInitialEquationID(int block_number, int equation_number) const
   {
-    return eq_idx_orig2block[equation_number] - get<1>(block_type_firstequation_size_mfs[block_number]);
+    return eq_idx_orig2block[equation_number] - getBlockFirstEquation(block_number);
   };
   //! Return the position of variable_number in the block number belonging to the block block_number
   int
   getBlockInitialVariableID(int block_number, int variable_number) const
   {
-    return endo_idx_orig2block[variable_number] - get<1>(block_type_firstequation_size_mfs[block_number]);
+    return endo_idx_orig2block[variable_number] - getBlockFirstEquation(block_number);
   };
   //! Return the position of variable_number in the block number belonging to the block block_number
   virtual int getBlockInitialExogenousID(int block_number, int variable_number) const = 0;
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index eceaa37d..1da00b22 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -703,7 +703,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
                                block_size,
                                endo_idx_block2orig,
                                eq_idx_block2orig,
-                               blocks_linear[block],
+                               blocks[block].linear,
                                symbol_table.endo_nbr(),
                                0,
                                0,
@@ -1142,9 +1142,9 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
 
       cout << "Finding the optimal block decomposition of the model ..." << endl;
 
-      auto [blocks, variable_lag_lead, n_static, n_forward, n_backward, n_mixed] = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, equation_type_and_normalized_equation, false);
+      auto variable_lag_lead = computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, false);
 
-      reduceBlocksAndTypeDetermination(blocks, equation_type_and_normalized_equation, n_static, n_forward, n_backward, n_mixed, false);
+      reduceBlocksAndTypeDetermination(false);
 
       printBlockDecomposition();