diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index a4ae1b2e021c82486543da05c01faf9cd21d0053..a877ce0b0ffd71abb45de1c85416ef3f3767a329 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -187,7 +187,7 @@ DynamicModel::computeTemporaryTermsOrdered()
   v_temporary_terms.clear();
   map_idx.clear();
 
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
   v_temporary_terms = vector<vector<temporary_terms_t>>(nb_blocks);
   v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
   temporary_terms.clear();
@@ -198,9 +198,8 @@ DynamicModel::computeTemporaryTermsOrdered()
         {
           reference_count.clear();
           temporary_terms.clear();
-          int block_size = getBlockSize(block);
-          int block_nb_mfs = getBlockMfs(block);
-          int block_nb_recursives = block_size - block_nb_mfs;
+          int block_size = blocks[block].size;
+          int block_nb_recursives = blocks[block].getRecursiveSize();
           v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
           for (int i = 0; i < block_size; i++)
             {
@@ -229,9 +228,8 @@ DynamicModel::computeTemporaryTermsOrdered()
       for (int block = 0; block < nb_blocks; block++)
         {
           // Compute the temporary terms reordered
-          int block_size = getBlockSize(block);
-          int block_nb_mfs = getBlockMfs(block);
-          int block_nb_recursives = block_size - block_nb_mfs;
+          int block_size = blocks[block].size;
+          int block_nb_recursives = blocks[block].getRecursiveSize();
           v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
           for (int i = 0; i < block_size; i++)
             {
@@ -256,9 +254,8 @@ DynamicModel::computeTemporaryTermsOrdered()
       for (int block = 0; block < nb_blocks; block++)
         {
           // Collect the temporary terms reordered
-          int block_size = getBlockSize(block);
-          int block_nb_mfs = getBlockMfs(block);
-          int block_nb_recursives = block_size - block_nb_mfs;
+          int block_size = blocks[block].size;
+          int block_nb_recursives = blocks[block].getRecursiveSize();
           set<int> temporary_terms_in_use;
           for (int i = 0; i < block_size; i++)
             {
@@ -320,7 +317,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
 
   //----------------------------------------------------------------------
   //For each block
-  for (int block = 0; block < getNbBlocks(); block++)
+  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
     {
 
       //recursive_variables.clear();
@@ -330,10 +327,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       nze_other_endo = derivative_other_endo[block].size();
       nze_exo = derivative_exo[block].size();
       nze_exo_det = derivative_exo_det[block].size();
-      BlockSimulationType simulation_type = getBlockSimulationType(block);
-      int block_size = getBlockSize(block);
-      int block_mfs = getBlockMfs(block);
-      int block_recursive = block_size - block_mfs;
+      BlockSimulationType simulation_type = blocks[block].simulation_type;
+      int block_size = blocks[block].size;
+      int block_mfs = blocks[block].mfs_size;
+      int block_recursive = blocks[block].getRecursiveSize();
       deriv_node_temp_terms_t tef_terms;
       local_output_type = ExprNodeOutputType::matlabDynamicModelSparse;
       if (global_temporary_terms)
@@ -439,13 +436,13 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                 || simulation_type == BlockSimulationType::solveBackwardSimple
                 || simulation_type == BlockSimulationType::evaluateBackward
                 || simulation_type == BlockSimulationType::evaluateForward)
-               && getBlockFirstEquation(block) < prologue)
+               && blocks[block].first_equation < prologue)
         block_type = BlockType::prologue;
       else if ((simulation_type == BlockSimulationType::solveForwardSimple
                 || simulation_type == BlockSimulationType::solveBackwardSimple
                 || simulation_type == BlockSimulationType::evaluateBackward
                 || simulation_type == BlockSimulationType::evaluateForward)
-               && getBlockFirstEquation(block) >= static_cast<int>(equations.size()) - epilogue)
+               && blocks[block].first_equation >= static_cast<int>(equations.size()) - epilogue)
         block_type = BlockType::epilogue;
       else
         block_type = BlockType::simultans;
@@ -1161,7 +1158,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
   FDIMT_ fdimt(temporary_terms.size());
   fdimt.write(code_file, instruction_number);
 
-  for (int block = 0; block < getNbBlocks(); block++)
+  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
     {
       feedback_variables.clear();
       if (block > 0)
@@ -1171,10 +1168,10 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
         }
       int count_u;
       int u_count_int = 0;
-      BlockSimulationType simulation_type = getBlockSimulationType(block);
-      int block_size = getBlockSize(block);
-      int block_mfs = getBlockMfs(block);
-      int block_recursive = block_size - block_mfs;
+      BlockSimulationType simulation_type = blocks[block].simulation_type;
+      int block_size = blocks[block].size;
+      int block_mfs = blocks[block].mfs_size;
+      int block_recursive = blocks[block].getRecursiveSize();
       int block_max_lag = max_leadlag_block[block].first;
       int block_max_lead = max_leadlag_block[block].second;
 
@@ -1245,7 +1242,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
 
       FBEGINBLOCK_ fbeginblock(block_mfs,
                                simulation_type,
-                               getBlockFirstEquation(block),
+                               blocks[block].first_equation,
                                block_size,
                                endo_idx_block2orig,
                                eq_idx_block2orig,
@@ -1808,9 +1805,9 @@ DynamicModel::Write_Inf_To_Bin_File_Block(const string &basename, int num,
       exit(EXIT_FAILURE);
     }
   u_count_int = 0;
-  int block_size = getBlockSize(num);
-  int block_mfs = getBlockMfs(num);
-  int block_recursive = block_size - block_mfs;
+  int block_size = blocks[num].size;
+  int block_mfs = blocks[num].mfs_size;
+  int block_recursive = blocks[num].getRecursiveSize();
   for (const auto &[eq, var, lag, ignore] : blocks_derivatives[num])
     {
       if (lag != 0 && !is_two_boundaries)
@@ -1911,14 +1908,13 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                     << "    ys=y(it_,:);" << endl;
   tmp.str("");
   tmp_eq.str("");
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
   int block = 0;
   for (int count_call = 1; block < nb_blocks; block++, count_call++)
     {
-      int block_size = getBlockSize(block),
-        block_mfs = getBlockMfs(block),
-        block_recursive = block_size - block_mfs;
-      BlockSimulationType simulation_type = getBlockSimulationType(block);
+      int block_size = blocks[block].size,
+        block_recursive = blocks[block].getRecursiveSize();
+      BlockSimulationType simulation_type = blocks[block].simulation_type;
 
       if (simulation_type == BlockSimulationType::evaluateForward
           || simulation_type == BlockSimulationType::evaluateBackward)
@@ -2006,10 +2002,9 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                     << "  oo_.deterministic_simulation.status = 0;" << endl;
   for (block = 0; block < nb_blocks; block++)
     {
-      int block_size = getBlockSize(block);
-      int block_mfs = getBlockMfs(block);
-      int block_recursive = block_size - block_mfs;
-      BlockSimulationType simulation_type = getBlockSimulationType(block);
+      int block_size = blocks[block].size;
+      int block_recursive = blocks[block].getRecursiveSize();
+      BlockSimulationType simulation_type = blocks[block].simulation_type;
 
       if (simulation_type == BlockSimulationType::evaluateForward && block_size)
         {
@@ -3135,13 +3130,13 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
       vector<int> state_equ;
       int count_lead_lag_incidence = 0;
       int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo, max_lag_exo_det, max_lead_exo_det;
-      int nb_blocks = getNbBlocks();
+      int nb_blocks = blocks.size();
       for (int block = 0; block < nb_blocks; block++)
         {
           //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
           count_lead_lag_incidence = 0;
-          BlockSimulationType simulation_type = getBlockSimulationType(block);
-          int block_size = getBlockSize(block);
+          BlockSimulationType simulation_type = blocks[block].simulation_type;
+          int block_size = blocks[block].size;
           max_lag = max_leadlag_block[block].first;
           max_lead = max_leadlag_block[block].second;
           max_lag_endo = endo_max_leadlag_block[block].first;
@@ -3180,7 +3175,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                  << "block_structure.block(" << block+1 << ").maximum_exo_det_lag = " << max_lag_exo_det << ";" << endl
                  << "block_structure.block(" << block+1 << ").maximum_exo_det_lead = " << max_lead_exo_det << ";" << endl
                  << "block_structure.block(" << block+1 << ").endo_nbr = " << block_size << ";" << endl
-                 << "block_structure.block(" << block+1 << ").mfs = " << getBlockMfs(block) << ";" << endl
+                 << "block_structure.block(" << block+1 << ").mfs = " << blocks[block].mfs_size << ";" << endl
                  << "block_structure.block(" << block+1 << ").equation = [" << tmp_s_eq.str() << "];" << endl
                  << "block_structure.block(" << block+1 << ").variable = [" << tmp_s.str() << "];" << endl
                  << "block_structure.block(" << block+1 << ").exo_nbr = " << getBlockExoSize(block) << ";" << endl
@@ -3220,7 +3215,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                 bool OK = true;
                 int j;
                 for (j = 0; j < block && OK; j++)
-                  for (int k = 0; k < getBlockSize(j) && OK; k++)
+                  for (int k = 0; k < blocks[j].size && OK; k++)
                     OK = other_endo != getBlockVariableID(j, k);
                 if (!OK)
                   output << " " << j;
@@ -3408,7 +3403,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
 
           for (int block = 0; block < nb_blocks; block++)
             {
-              int block_size = getBlockSize(block);
+              int block_size = blocks[block].size;
               int nze = 0;
 
               for (int i = 0; i < block_size; i++)
@@ -4857,9 +4852,9 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
       equation_block.resize(equations.size());
       variable_block_lead_lag.clear();
       variable_block_lead_lag.resize(equations.size());
-      for (int i = 0; i < getNbBlocks(); i++)
+      for (int i = 0; i < static_cast<int>(blocks.size()); i++)
         {
-          for (int j = 0; j < getBlockSize(i); j++)
+          for (int j = 0; j < blocks[i].size; j++)
             {
               equation_block[eq_idx_block2orig[k]] = i;
               int l = endo_idx_block2orig[k];
@@ -4983,21 +4978,22 @@ DynamicModel::get_Derivatives(int block)
 {
   int max_lag, max_lead;
   map<tuple<int, int, int, int, int>, int> Derivatives;
-  BlockSimulationType simulation_type = getBlockSimulationType(block);
+  BlockSimulationType simulation_type = blocks[block].simulation_type;
   if (simulation_type == BlockSimulationType::evaluateBackward
       || simulation_type == BlockSimulationType::evaluateForward)
     {
       max_lag = 1;
       max_lead = 1;
-      setBlockLeadLag(block, max_lag, max_lead);
+      blocks[block].max_lag = max_lag;
+      blocks[block].max_lead = max_lead;
     }
   else
     {
-      max_lag = getBlockMaxLag(block);
-      max_lead = getBlockMaxLead(block);
+      max_lag = blocks[block].max_lag;
+      max_lead = blocks[block].max_lead;
     }
-  int block_size = getBlockSize(block);
-  int block_nb_recursive = block_size - getBlockMfs(block);
+  int block_size = blocks[block].size;
+  int block_nb_recursive = blocks[block].getRecursiveSize();
   for (int lag = -max_lag; lag <= max_lead; lag++)
     {
       for (int eq = 0; eq < block_size; eq++)
@@ -5045,15 +5041,13 @@ void
 DynamicModel::computeChainRuleJacobian()
 {
   map<int, expr_t> recursive_variables;
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
   blocks_derivatives.resize(nb_blocks);
   for (int block = 0; block < nb_blocks; block++)
     {
       block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
       recursive_variables.clear();
-      int block_size = getBlockSize(block);
-      int block_nb_mfs = getBlockMfs(block);
-      int block_nb_recursives = block_size - block_nb_mfs;
+      int block_nb_recursives = blocks[block].getRecursiveSize();
       for (int i = 0; i < block_nb_recursives; i++)
         {
           if (getBlockEquationType(block, i) == EquationType::evaluate_s)
@@ -5089,16 +5083,13 @@ DynamicModel::collect_block_first_order_derivatives()
 {
   //! vector for an equation or a variable indicates the block number
   vector<int> equation_2_block(eq_idx_block2orig.size()), variable_2_block(endo_idx_block2orig.size());
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
   for (int block = 0; block < nb_blocks; block++)
-    {
-      int block_size = getBlockSize(block);
-      for (int i = 0; i < block_size; i++)
-        {
-          equation_2_block[getBlockEquationID(block, i)] = block;
-          variable_2_block[getBlockVariableID(block, i)] = block;
-        }
-    }
+    for (int i = 0; i < blocks[block].size; i++)
+      {
+        equation_2_block[getBlockEquationID(block, i)] = block;
+        variable_2_block[getBlockVariableID(block, i)] = block;
+      }
   other_endo_block = vector<lag_var_t>(nb_blocks);
   exo_block = vector<lag_var_t>(nb_blocks);
   exo_det_block = vector<lag_var_t>(nb_blocks);
@@ -5227,7 +5218,7 @@ DynamicModel::collect_block_first_order_derivatives()
 void
 DynamicModel::collectBlockVariables()
 {
-  for (int block = 0; block < getNbBlocks(); block++)
+  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
     {
       int prev_var = -1;
       int prev_lag = -999999999;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 465d2f525da01644e4bfc66fc38dbe0c1eb2cc4e..b4bad7ff46a63904d0656a6737f8becbad618caf 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -882,19 +882,19 @@ void
 ModelTree::printBlockDecomposition() const
 {
   int largest_block = 0, Nb_SimulBlocks = 0, Nb_feedback_variable = 0;
-  int Nb_TotalBlocks = getNbBlocks();
+  int Nb_TotalBlocks = blocks.size();
   for (int block = 0; block < Nb_TotalBlocks; block++)
-    if (BlockSimulationType simulation_type = getBlockSimulationType(block);
+    if (BlockSimulationType simulation_type = blocks[block].simulation_type;
         simulation_type == BlockSimulationType::solveForwardComplete
         || simulation_type == BlockSimulationType::solveBackwardComplete
         || simulation_type == BlockSimulationType::solveTwoBoundariesComplete)
       {
         Nb_SimulBlocks++;
-        if (int size = getBlockSize(block);
+        if (int size = blocks[block].size;
             size > largest_block)
           {
             largest_block = size;
-            Nb_feedback_variable = getBlockMfs(block);
+            Nb_feedback_variable = blocks[block].mfs_size;
           }
       }
 
@@ -1051,12 +1051,12 @@ void
 ModelTree::determineLinearBlocks()
 {
   // Note that field “linear” in class BlockInfo defaults to true
-  for (int block = 0; block < getNbBlocks(); block++)
+  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
     {
-      BlockSimulationType simulation_type = getBlockSimulationType(block);
-      int block_size = getBlockSize(block);
+      BlockSimulationType simulation_type = blocks[block].simulation_type;
+      int block_size = blocks[block].size;
       block_derivatives_equation_variable_laglead_nodeid_t derivatives_block = blocks_derivatives[block];
-      int first_variable_position = getBlockFirstEquation(block);
+      int first_variable_position = blocks[block].first_equation;
       if (simulation_type == BlockSimulationType::solveBackwardComplete
           || simulation_type == BlockSimulationType::solveForwardComplete)
         for (const auto &[ignore, ignore2, lag, d1] : derivatives_block)
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 17ad90a7a49d5712fb3110030a52cb2347e01259..34862c3cbc946acd7d705b79faa4e0fd13145407 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -208,6 +208,8 @@ protected:
     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};
+
+    inline int getRecursiveSize() const { return size - mfs_size; };
   };
 
   // Stores various informations on the blocks
@@ -337,108 +339,59 @@ protected:
                                       vector<BinaryOpNode *> &equations, vector<int> &equations_lineno,
                                       EquationTags &equation_tags, bool static_equations) const;
 
-  //! Return the number of blocks
-  int
-  getNbBlocks() const
-  {
-    return blocks.size();
-  };
-  //! Determine the simulation type of each block
-  BlockSimulationType
-  getBlockSimulationType(int block_number) const
-  {
-    return blocks[block_number].simulation_type;
-  };
-  //! Return the first equation number of a block
-  int
-  getBlockFirstEquation(int block_number) const
-  {
-    return blocks[block_number].first_equation;
-  };
-  //! Return the size of the block block_number
-  int
-  getBlockSize(int block_number) const
-  {
-    return blocks[block_number].size;
-  };
   //! Return the number of exogenous variable in the block block_number
   virtual int getBlockExoSize(int block_number) const = 0;
   //! Return the number of colums in the jacobian matrix for exogenous variable in the block block_number
   virtual int getBlockExoColSize(int block_number) const = 0;
-  //! Return the number of feedback variable of the block block_number
-  int
-  getBlockMfs(int block_number) const
-  {
-    return blocks[block_number].mfs_size;
-  };
-  //! Return the maximum lag in a block
-  int
-  getBlockMaxLag(int block_number) const
-  {
-    return blocks[block_number].max_lag;
-  };
-  //! Return the maximum lead in a block
-  int
-  getBlockMaxLead(int block_number) const
-  {
-    return blocks[block_number].max_lead;
-  };
-  inline void
-  setBlockLeadLag(int block, int max_lag, int 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
+  //! Return the type of equation belonging to the block
   EquationType
-  getBlockEquationType(int block_number, int equation_number) const
+  getBlockEquationType(int blk, int eq) const
   {
-    return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].first;
+    return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation+eq]].first;
   };
   //! Return true if the equation has been normalized
   bool
-  isBlockEquationRenormalized(int block_number, int equation_number) const
+  isBlockEquationRenormalized(int blk, int eq) const
   {
-    return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].first == EquationType::evaluate_s;
+    return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]].first == EquationType::evaluate_s;
   };
-  //! Return the expr_t of the equation equation_number belonging to the block block_number
+  //! Return the expr_t of equation belonging to the block
   expr_t
-  getBlockEquationExpr(int block_number, int equation_number) const
+  getBlockEquationExpr(int blk, int eq) const
   {
-    return equations[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]];
+    return equations[eq_idx_block2orig[blocks[blk].first_equation + eq]];
   };
-  //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
+  //! Return the expr_t of renormalized equation belonging to the block
   expr_t
-  getBlockEquationRenormalizedExpr(int block_number, int equation_number) const
+  getBlockEquationRenormalizedExpr(int blk, int eq) const
   {
-    return equation_type_and_normalized_equation[eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number]].second;
+    return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]].second;
   };
-  //! Return the original number of equation equation_number belonging to the block block_number
+  //! Return the original number of equation belonging to the block
   int
-  getBlockEquationID(int block_number, int equation_number) const
+  getBlockEquationID(int blk, int eq) const
   {
-    return eq_idx_block2orig[getBlockFirstEquation(block_number)+equation_number];
+    return eq_idx_block2orig[blocks[blk].first_equation + eq];
   };
-  //! Return the original number of variable variable_number belonging to the block block_number
+  //! Return the original number of variable belonging to the block
   int
-  getBlockVariableID(int block_number, int variable_number) const
+  getBlockVariableID(int blk, int var) const
   {
-    return endo_idx_block2orig[getBlockFirstEquation(block_number)+variable_number];
+    return endo_idx_block2orig[blocks[blk].first_equation + var];
   };
   //! 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;
-  //! Return the position of equation_number in the block number belonging to the block block_number
+  //! Return the position of an equation (given by its original index) inside its block
   int
-  getBlockInitialEquationID(int block_number, int equation_number) const
+  getBlockInitialEquationID(int blk, int eq) const
   {
-    return eq_idx_orig2block[equation_number] - getBlockFirstEquation(block_number);
+    return eq_idx_orig2block[eq] - blocks[blk].first_equation;
   };
-  //! Return the position of variable_number in the block number belonging to the block block_number
+  //! Return the position of a variable (given by its original index) inside its block
   int
-  getBlockInitialVariableID(int block_number, int variable_number) const
+  getBlockInitialVariableID(int blk, int var) const
   {
-    return endo_idx_orig2block[variable_number] - getBlockFirstEquation(block_number);
+    return endo_idx_orig2block[var] - blocks[blk].first_equation;
   };
   //! 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 04ff1aeb0947a9b724717dce1c0889b2c2170d5f..503ce8f33a6db0924f50e091516142a936bc234b 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -161,7 +161,7 @@ StaticModel::computeTemporaryTermsOrdered()
   ostringstream tmp_s;
   map_idx.clear();
 
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
   v_temporary_terms = vector< vector<temporary_terms_t>>(nb_blocks);
   v_temporary_terms_local = vector< vector<temporary_terms_t>>(nb_blocks);
 
@@ -178,9 +178,8 @@ StaticModel::computeTemporaryTermsOrdered()
       map<expr_t, pair<int, int>> first_occurence_local;
       temporary_terms_t temporary_terms_l;
 
-      int block_size = getBlockSize(block);
-      int block_nb_mfs = getBlockMfs(block);
-      int block_nb_recursives = block_size - block_nb_mfs;
+      int block_size = blocks[block].size;
+      int block_nb_recursives = blocks[block].getRecursiveSize();
       v_temporary_terms_local[block] = vector<temporary_terms_t>(block_size);
 
       for (int i = 0; i < block_size; i++)
@@ -206,9 +205,8 @@ StaticModel::computeTemporaryTermsOrdered()
   for (int block = 0; block < nb_blocks; block++)
     {
       // Compute the temporary terms reordered
-      int block_size = getBlockSize(block);
-      int block_nb_mfs = getBlockMfs(block);
-      int block_nb_recursives = block_size - block_nb_mfs;
+      int block_size = blocks[block].size;
+      int block_nb_recursives = blocks[block].getRecursiveSize();
       v_temporary_terms[block] = vector<temporary_terms_t>(block_size);
       for (int i = 0; i < block_size; i++)
         {
@@ -230,9 +228,8 @@ StaticModel::computeTemporaryTermsOrdered()
   for (int block = 0; block < nb_blocks; block++)
     {
       // Collecte the temporary terms reordered
-      int block_size = getBlockSize(block);
-      int block_nb_mfs = getBlockMfs(block);
-      int block_nb_recursives = block_size - block_nb_mfs;
+      int block_size = blocks[block].size;
+      int block_nb_recursives = blocks[block].getRecursiveSize();
       set<int> temporary_terms_in_use;
       for (int i = 0; i < block_size; i++)
         {
@@ -249,7 +246,7 @@ StaticModel::computeTemporaryTermsOrdered()
           expr_t id = get<3>(it);
           id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
         }
-      for (int i = 0; i < getBlockSize(block); i++)
+      for (int i = 0; i < blocks[block].size; i++)
         for (const auto &it : v_temporary_terms[block][i])
           it->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
       v_temporary_terms_inuse[block] = temporary_terms_in_use;
@@ -286,15 +283,15 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
 
   //----------------------------------------------------------------------
   //For each block
-  for (int block = 0; block < getNbBlocks(); block++)
+  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
     {
       //recursive_variables.clear();
       feedback_variables.clear();
       //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
-      BlockSimulationType simulation_type = getBlockSimulationType(block);
-      int block_size = getBlockSize(block);
-      int block_mfs = getBlockMfs(block);
-      int block_recursive = block_size - block_mfs;
+      BlockSimulationType simulation_type = blocks[block].simulation_type;
+      int block_size = blocks[block].size;
+      int block_mfs = blocks[block].mfs_size;
+      int block_recursive = blocks[block].getRecursiveSize();
 
       tmp1_output.str("");
       tmp1_output << packageDir(basename + ".block") << "/static_" << block+1 << ".m";
@@ -319,13 +316,13 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
                 || simulation_type == BlockSimulationType::solveBackwardSimple
                 || simulation_type == BlockSimulationType::evaluateBackward
                 || simulation_type == BlockSimulationType::evaluateForward)
-               && getBlockFirstEquation(block) < prologue)
+               && blocks[block].first_equation < prologue)
         block_type = BlockType::prologue;
       else if ((simulation_type == BlockSimulationType::solveForwardSimple
                 || simulation_type == BlockSimulationType::solveBackwardSimple
                 || simulation_type == BlockSimulationType::evaluateBackward
                 || simulation_type == BlockSimulationType::evaluateForward)
-               && getBlockFirstEquation(block) >= static_cast<int>(equations.size()) - epilogue)
+               && blocks[block].first_equation >= static_cast<int>(equations.size()) - epilogue)
         block_type = BlockType::epilogue;
       else
         block_type = BlockType::simultans;
@@ -673,7 +670,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
   FDIMST_ fdimst(temporary_terms.size());
   fdimst.write(code_file, instruction_number);
 
-  for (int block = 0; block < getNbBlocks(); block++)
+  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
     {
       feedback_variables.clear();
       if (block > 0)
@@ -683,10 +680,10 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
         }
       int count_u;
       int u_count_int = 0;
-      BlockSimulationType simulation_type = getBlockSimulationType(block);
-      int block_size = getBlockSize(block);
-      int block_mfs = getBlockMfs(block);
-      int block_recursive = block_size - block_mfs;
+      BlockSimulationType simulation_type = blocks[block].simulation_type;
+      int block_size = blocks[block].size;
+      int block_mfs = blocks[block].mfs_size;
+      int block_recursive = blocks[block].getRecursiveSize();
 
       if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple
           || simulation_type == BlockSimulationType::solveTwoBoundariesComplete
@@ -699,7 +696,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
 
       FBEGINBLOCK_ fbeginblock(block_mfs,
                                simulation_type,
-                               getBlockFirstEquation(block),
+                               blocks[block].first_equation,
                                block_size,
                                endo_idx_block2orig,
                                eq_idx_block2orig,
@@ -1055,9 +1052,9 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &basename, int num,
       exit(EXIT_FAILURE);
     }
   u_count_int = 0;
-  int block_size = getBlockSize(num);
-  int block_mfs = getBlockMfs(num);
-  int block_recursive = block_size - block_mfs;
+  int block_size = blocks[num].size;
+  int block_mfs = blocks[num].mfs_size;
+  int block_recursive = blocks[num].getRecursiveSize();
   for (const auto &[eq, var, ignore, ignore2] : blocks_derivatives[num])
     {
       int lag = 0;
@@ -1970,7 +1967,7 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const
          << "  var_index = [];" << endl << endl
          << "  switch nblock" << endl;
 
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
 
   for (int b = 0; b < nb_blocks; b++)
     {
@@ -1978,14 +1975,14 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const
 
       output << "    case " << b+1 << endl;
 
-      BlockSimulationType simulation_type = getBlockSimulationType(b);
+      BlockSimulationType simulation_type = blocks[b].simulation_type;
 
       if (simulation_type == BlockSimulationType::evaluateBackward
           || simulation_type == BlockSimulationType::evaluateForward)
         {
           output << "      y_tmp = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl;
           ostringstream tmp;
-          for (int i = 0; i < getBlockSize(b); i++)
+          for (int i = 0; i < blocks[b].size; i++)
             tmp << " " << getBlockVariableID(b, i)+1;
           output << "      var_index = [" << tmp.str() << "];" << endl
                  << "      residual  = y(var_index) - y_tmp(var_index);" << endl
@@ -2011,11 +2008,11 @@ StaticModel::writeOutput(ostream &output, bool block) const
   if (!block)
     return;
 
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
   for (int b = 0; b < nb_blocks; b++)
     {
-      BlockSimulationType simulation_type = getBlockSimulationType(b);
-      int block_size = getBlockSize(b);
+      BlockSimulationType simulation_type = blocks[b].simulation_type;
+      int block_size = blocks[b].size;
       ostringstream tmp_s, tmp_s_eq;
       tmp_s.str("");
       tmp_s_eq.str("");
@@ -2026,7 +2023,7 @@ StaticModel::writeOutput(ostream &output, bool block) const
         }
       output << "block_structure_stat.block(" << b+1 << ").Simulation_Type = " << static_cast<int>(simulation_type) << ";" << endl
              << "block_structure_stat.block(" << b+1 << ").endo_nbr = " << block_size << ";" << endl
-             << "block_structure_stat.block(" << b+1 << ").mfs = " << getBlockMfs(b) << ";" << endl
+             << "block_structure_stat.block(" << b+1 << ").mfs = " << blocks[b].mfs_size << ";" << endl
              << "block_structure_stat.block(" << b+1 << ").equation = [" << tmp_s_eq.str() << "];" << endl
              << "block_structure_stat.block(" << b+1 << ").variable = [" << tmp_s.str() << "];" << endl;
     }
@@ -2111,8 +2108,8 @@ map<tuple<int, int, int, int, int>, int>
 StaticModel::get_Derivatives(int block)
 {
   map<tuple<int, int, int, int, int>, int> Derivatives;
-  int block_size = getBlockSize(block);
-  int block_nb_recursive = block_size - getBlockMfs(block);
+  int block_size = blocks[block].size;
+  int block_nb_recursive = blocks[block].getRecursiveSize();
   int lag = 0;
   for (int eq = 0; eq < block_size; eq++)
     {
@@ -2159,16 +2156,15 @@ void
 StaticModel::computeChainRuleJacobian()
 {
   map<int, expr_t> recursive_variables;
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
   blocks_derivatives.resize(nb_blocks);
   for (int block = 0; block < nb_blocks; block++)
     {
       block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
       recursive_variables.clear();
-      BlockSimulationType simulation_type = getBlockSimulationType(block);
-      int block_size = getBlockSize(block);
-      int block_nb_mfs = getBlockMfs(block);
-      int block_nb_recursives = block_size - block_nb_mfs;
+      BlockSimulationType simulation_type = blocks[block].simulation_type;
+      int block_size = blocks[block].size;
+      int block_nb_recursives = blocks[block].getRecursiveSize();
       if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
           || simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
         {
@@ -2231,10 +2227,10 @@ StaticModel::collect_block_first_order_derivatives()
 {
   //! vector for an equation or a variable indicates the block number
   vector<int> equation_2_block(eq_idx_block2orig.size()), variable_2_block(endo_idx_block2orig.size());
-  int nb_blocks = getNbBlocks();
+  int nb_blocks = blocks.size();
   for (int block = 0; block < nb_blocks; block++)
     {
-      int block_size = getBlockSize(block);
+      int block_size = blocks[block].size;
       for (int i = 0; i < block_size; i++)
         {
           equation_2_block[getBlockEquationID(block, i)] = block;