diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 6866a8128f71816568ea53c6d293d9d8ca779a46..4db74ef2e5a588553ec0d1dcf2e123b03dbc0513 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -106,7 +106,6 @@ DynamicModel::DynamicModel(const DynamicModel &m) :
   blocks_jacob_cols_other_endo{m.blocks_jacob_cols_other_endo},
   blocks_jacob_cols_exo{m.blocks_jacob_cols_exo},
   blocks_jacob_cols_exo_det{m.blocks_jacob_cols_exo_det},
-  global_temporary_terms{m.global_temporary_terms},
   var_expectation_functions_to_write{m.var_expectation_functions_to_write},
   pac_mce_alpha_symb_ids{m.pac_mce_alpha_symb_ids},
   pac_h0_indices{m.pac_h0_indices},
@@ -168,7 +167,6 @@ DynamicModel::operator=(const DynamicModel &m)
   blocks_jacob_cols_other_endo = m.blocks_jacob_cols_other_endo;
   blocks_jacob_cols_exo = m.blocks_jacob_cols_exo;
   blocks_jacob_cols_exo_det = m.blocks_jacob_cols_exo_det;
-  global_temporary_terms = m.global_temporary_terms;
 
   var_expectation_functions_to_write = m.var_expectation_functions_to_write;
 
@@ -187,11 +185,11 @@ DynamicModel::operator=(const DynamicModel &m)
 }
 
 void
-DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const
+DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag) const
 {
   if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) });
       it != derivatives[1].end())
-    it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+    it->second->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
   else
     {
       FLDZ_ fldz;
@@ -200,11 +198,11 @@ DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_n
 }
 
 void
-DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const map_idx_t &map_idx) const
+DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag) const
 {
   if (auto it = blocks_derivatives[blk].find({ eq, var, lag });
       it != blocks_derivatives[blk].end())
-    it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+    it->second->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
   else
     {
       FLDZ_ fldz;
@@ -213,106 +211,16 @@ DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &inst
 }
 
 void
-DynamicModel::computeTemporaryTermsOrdered()
+DynamicModel::additionalBlockTemporaryTerms(int blk,
+                                            vector<temporary_terms_t> &blocks_temporary_terms,
+                                            map<expr_t, pair<int, int>> &reference_count) const
 {
-  map<expr_t, pair<int, int>> first_occurence;
-  map<expr_t, int> reference_count;
-  BinaryOpNode *eq_node;
-  ostringstream tmp_s;
-  v_temporary_terms.clear();
-  map_idx.clear();
-
-  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();
-
-  if (!global_temporary_terms)
-    {
-      for (int block = 0; block < nb_blocks; block++)
-        {
-          reference_count.clear();
-          temporary_terms.clear();
-          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++)
-            {
-              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
-              else
-                {
-                  eq_node = getBlockEquationExpr(block, i);
-                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
-                }
-            }
-          for (const auto &[ignore, d] : blocks_derivatives[block])
-            d->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
-          for (const auto &[ignore, d] : blocks_derivatives_other_endo[block])
-            d->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
-          v_temporary_terms_inuse[block] = {};
-        }
-    }
-  else
-    {
-      for (int block = 0; block < nb_blocks; block++)
-        {
-          // Compute the temporary terms reordered
-          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++)
-            {
-              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
-              else
-                {
-                  eq_node = getBlockEquationExpr(block, i);
-                  eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
-                }
-            }
-          for (const auto &[ignore, d] : blocks_derivatives[block])
-            d->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
-          for (const auto &[ignore, d] : blocks_derivatives_other_endo[block])
-            d->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
-        }
-      for (int block = 0; block < nb_blocks; block++)
-        {
-          // Collect the temporary terms reordered
-          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++)
-            {
-              if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-              else
-                {
-                  eq_node = getBlockEquationExpr(block, i);
-                  eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-                }
-            }
-          for (const auto &[ignore, d] : blocks_derivatives[block])
-            d->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          for (const auto &[ignore, d] : blocks_derivatives_other_endo[block])
-            d->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          for (const auto &[ignore, d] : blocks_derivatives_exo[block])
-            d->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          for (const auto &[ignore, d] : blocks_derivatives_exo_det[block])
-            d->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          v_temporary_terms_inuse[block] = temporary_terms_in_use;
-        }
-      computeTemporaryTermsMapping();
-    }
-}
-
-void
-DynamicModel::computeTemporaryTermsMapping()
-{
-  // Add a mapping form node ID to temporary terms order
-  int j = 0;
-  for (auto temporary_term : temporary_terms)
-    map_idx[temporary_term->idx] = j++;
+  for (const auto &[ignore, d] : blocks_derivatives_exo[blk])
+    d->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+  for (const auto &[ignore, d] : blocks_derivatives_exo_det[blk])
+    d->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+  for (const auto &[ignore, d] : blocks_derivatives_other_endo[blk])
+    d->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
 }
 
 void
@@ -325,15 +233,12 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
   ostringstream Ufoss;
   vector<string> Uf(symbol_table.endo_nbr(), "");
   map<expr_t, int> reference_count;
-  temporary_terms_t local_temporary_terms;
   ofstream output;
   int nze, nze_exo, nze_exo_det, nze_other_endo;
   vector<int> feedback_variables;
-  ExprNodeOutputType local_output_type;
 
-  local_output_type = ExprNodeOutputType::matlabDynamicModelSparse;
-  if (global_temporary_terms)
-    local_temporary_terms = temporary_terms;
+  temporary_terms_t temporary_terms; // Temp terms written so far
+  ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabDynamicModelSparse;
 
   //----------------------------------------------------------------------
   //For each block
@@ -351,10 +256,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       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)
-        local_temporary_terms = temporary_terms;
 
       tmp1_output.str("");
       tmp1_output << packageDir(basename + ".block") << "/dynamic_" << block+1 << ".m";
@@ -416,31 +318,36 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
         }
 
       output << "  g2=0;g3=0;" << endl;
-      if (v_temporary_terms_inuse[block].size())
-        {
-          tmp_output.str("");
-          for (int it : v_temporary_terms_inuse[block])
-            tmp_output << " T" << it;
-          output << "  global" << tmp_output.str() << ";" << endl;
-        }
+
+      // Declare global temp terms from this block and the previous ones
+      bool global_keyword_written = false;
+      for (int blk2 = 0; blk2 <= block; blk2++)
+        for (auto tt : blocks_temporary_terms[blk2])
+          {
+            if (!global_keyword_written)
+              {
+                output << "  global";
+                global_keyword_written = true;
+              }
+            output << " ";
+            // In the following, "Static" is used to avoid getting the "(it_)" subscripting
+            tt->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, blocks_temporary_terms[blk2], {});
+          }
+      if (global_keyword_written)
+        output << ";" << endl;
+
+      // Initialize temp terms of this block
       if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
           || simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
         {
-          temporary_terms_t tt2;
-          for (int i = 0; i < block_size; i++)
+          output << "  " << "% //Temporary variables initialization" << endl
+                 << "  " << "T_zeros = zeros(y_kmin+periods, 1);" << endl;
+          for (auto it : blocks_temporary_terms[block])
             {
-              if (v_temporary_terms[block][i].size() && global_temporary_terms)
-                {
-                  output << "  " << "% //Temporary variables initialization" << endl
-                         << "  " << "T_zeros = zeros(y_kmin+periods, 1);" << endl;
-                  for (auto it : v_temporary_terms[block][i])
-                    {
-                      output << "  ";
-                      // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-                      it->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, local_temporary_terms, {});
-                      output << " = T_zeros;" << endl;
-                    }
-                }
+              output << "  ";
+              // In the following, "Static" is used to avoid getting the "(it_)" subscripting
+              it->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, blocks_temporary_terms[block], {});
+              output << " = T_zeros;" << endl;
             }
         }
       if (simulation_type == BlockSimulationType::solveBackwardSimple
@@ -472,29 +379,25 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
           sps = "  ";
         else
           sps = "";
+
+      deriv_node_temp_terms_t tef_terms;
+      for (auto it : blocks_temporary_terms[block])
+        {
+          if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+            it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+
+          output << "  " <<  sps;
+          it->writeOutput(output, local_output_type, blocks_temporary_terms[block], {}, tef_terms);
+          output << " = ";
+          it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms);
+          temporary_terms.insert(it);
+          output << ";" << endl;
+        }
+
       // The equations
       temporary_terms_idxs_t temporary_terms_idxs;
       for (int i = 0; i < block_size; i++)
         {
-          temporary_terms_t tt2;
-          if (v_temporary_terms[block].size())
-            {
-              output << "  " << "% //Temporary variables" << endl;
-              for (auto it : v_temporary_terms[block][i])
-                {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
-                    it->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
-
-                  output << "  " <<  sps;
-                  it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
-                  output << " = ";
-                  it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
-                  // Insert current node into tt2
-                  tt2.insert(it);
-                  output << ";" << endl;
-                }
-            }
-
           int variable_ID = getBlockVariableID(block, i);
           int equation_ID = getBlockEquationID(block, i);
           EquationType equ_type = getBlockEquationType(block, i);
@@ -503,7 +406,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
           lhs = eq_node->arg1;
           rhs = eq_node->arg2;
           tmp_output.str("");
-          lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
+          lhs->writeOutput(tmp_output, local_output_type, temporary_terms, {});
           switch (simulation_type)
             {
             case BlockSimulationType::evaluateBackward:
@@ -517,7 +420,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                 {
                   output << tmp_output.str();
                   output << " = ";
-                  rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+                  rhs->writeOutput(output, local_output_type, temporary_terms, {});
                 }
               else if (equ_type == EquationType::evaluate_s)
                 {
@@ -525,15 +428,15 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                   output << " = ";
                   if (isBlockEquationRenormalized(block, i))
                     {
-                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+                      rhs->writeOutput(output, local_output_type, temporary_terms, {});
                       output << endl << "    ";
                       tmp_output.str("");
                       eq_node = getBlockEquationRenormalizedExpr(block, i);
                       lhs = eq_node->arg1;
                       rhs = eq_node->arg2;
-                      lhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+                      lhs->writeOutput(output, local_output_type, temporary_terms, {});
                       output << " = ";
-                      rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+                      rhs->writeOutput(output, local_output_type, temporary_terms, {});
                     }
                 }
               else
@@ -570,7 +473,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
             end:
               output << tmp_output.str();
               output << ") - (";
-              rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+              rhs->writeOutput(output, local_output_type, temporary_terms, {});
               output << ");" << endl;
 #ifdef CONDITION
               if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
@@ -596,7 +499,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
           auto [eq, var, lag] = indices;
           int jacob_col = blocks_jacob_cols_endo[block].at({ var, lag });
           output << "      g1(" << eq+1 << ", " << jacob_col+1 << ") = ";
-          d->writeOutput(output, local_output_type, local_temporary_terms, {});
+          d->writeOutput(output, local_output_type, temporary_terms, {});
           output << ";" << endl;
         }
       for (const auto &[indices, d] : blocks_derivatives_exo[block])
@@ -604,7 +507,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
           auto [eq, var, lag] = indices;
           int jacob_col = blocks_jacob_cols_exo[block].at({ var, lag });
           output << "      g1_x(" << eq+1 << ", " << jacob_col+1 << ") = ";
-          d->writeOutput(output, local_output_type, local_temporary_terms, {});
+          d->writeOutput(output, local_output_type, temporary_terms, {});
           output << ";" << endl;
         }
       for (const auto &[indices, d] : blocks_derivatives_exo_det[block])
@@ -612,7 +515,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
           auto [eq, var, lag] = indices;
           int jacob_col = blocks_jacob_cols_exo_det[block].at({ var, lag });
           output << "      g1_xd(" << eq+1 << ", " << jacob_col+1 << ") = ";
-          d->writeOutput(output, local_output_type, local_temporary_terms, {});
+          d->writeOutput(output, local_output_type, temporary_terms, {});
           output << ";" << endl;
         }
       for (const auto &[indices, d] : blocks_derivatives_other_endo[block])
@@ -620,7 +523,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
           auto [eq, var, lag] = indices;
           int jacob_col = blocks_jacob_cols_other_endo[block].at({ var, lag });
           output << "      g1_o(" << eq+1 << ", " << jacob_col+1 << ") = ";
-          d->writeOutput(output, local_output_type, local_temporary_terms, {});
+          d->writeOutput(output, local_output_type, temporary_terms, {});
           output << ";" << endl;
         }
       output << "      varargout{1}=g1_x;" << endl
@@ -647,7 +550,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               if (lag == 0)
                 {
                   output << "    g1(" << eq+1 << ", " << var+1-block_recursive << ") = ";
-                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
+                  id->writeOutput(output, local_output_type, temporary_terms, {});
                   output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
                          << "(" << lag
                          << ") " << varr+1
@@ -700,7 +603,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                     tmp_output << "     g1(" << eq+1-block_recursive << "+Per_J_, "
                                << var+1-block_recursive << "+y_size*(it_" << lag-1 << ")) = ";
                   output << " " << tmp_output.str();
-                  id->writeOutput(output, local_output_type, local_temporary_terms, {});
+                  id->writeOutput(output, local_output_type, temporary_terms, {});
                   output << ";";
                   output << " %2 variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
                          << "(" << lag << ") " << varr+1
@@ -749,7 +652,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
 }
 
 void
-DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &map_idx) const
+DynamicModel::writeModelEquationsCode(const string &basename) const
 {
 
   ostringstream tmp_output;
@@ -780,6 +683,11 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
   Write_Inf_To_Bin_File(basename + "/model/bytecode/dynamic.bin", u_count_int, file_open, simulation_type == BlockSimulationType::solveTwoBoundariesComplete, symbol_table.endo_nbr());
   file_open = true;
 
+  // Compute the union of temporary terms from residuals and 1st derivatives
+  temporary_terms_t temporary_terms = temporary_terms_derivatives[0];
+  copy(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end(),
+       inserter(temporary_terms, temporary_terms.end()));
+
   //Temporary variables declaration
   FDIMT_ fdimt(temporary_terms.size());
   fdimt.write(code_file, instruction_number);
@@ -865,9 +773,9 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
                            other_endo);
   fbeginblock.write(code_file, instruction_number);
 
-  compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
+  compileTemporaryTerms(code_file, instruction_number, true, false);
 
-  compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, true, false);
+  compileModelEquations(code_file, instruction_number, true, false);
 
   FENDEQU_ fendequ;
   fendequ.write(code_file, instruction_number);
@@ -894,7 +802,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
           if (!my_derivatives[eq].size())
             my_derivatives[eq].clear();
           my_derivatives[eq].emplace_back(var, lag, count_u);
-          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+          d1->compile(code_file, instruction_number, false, temporary_terms_idxs, true, false);
 
           FSTPU_ fstpu(count_u);
           fstpu.write(code_file, instruction_number);
@@ -956,7 +864,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
           prev_lag = lag;
           count_col_endo++;
         }
-      d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+      d1->compile(code_file, instruction_number, false, temporary_terms_idxs, true, false);
       FSTPG3_ fstpg3(eq, var, lag, count_col_endo-1);
       fstpg3.write(code_file, instruction_number);
     }
@@ -975,7 +883,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
           prev_lag = lag;
           count_col_exo++;
         }
-      d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+      d1->compile(code_file, instruction_number, false, temporary_terms_idxs, true, false);
       FSTPG3_ fstpg3(eq, var, lag, count_col_exo-1);
       fstpg3.write(code_file, instruction_number);
     }
@@ -994,7 +902,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
 }
 
 void
-DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx, bool linear_decomposition) const
+DynamicModel::writeModelEquationsCode_Block(const string &basename, bool linear_decomposition) const
 {
   struct Uff_l
   {
@@ -1016,7 +924,6 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
   BinaryOpNode *eq_node;
   Uff Uf[symbol_table.endo_nbr()];
   map<expr_t, int> reference_count;
-  deriv_node_temp_terms_t tef_terms;
   vector<int> feedback_variables;
   bool file_open = false;
   string main_name;
@@ -1033,7 +940,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
     }
   //Temporary variables declaration
 
-  FDIMT_ fdimt(temporary_terms.size());
+  FDIMT_ fdimt(blocks_temporary_terms_idxs.size());
   fdimt.write(code_file, instruction_number);
 
   for (int block = 0; block < static_cast<int>(blocks.size()); block++)
@@ -1087,45 +994,33 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
       fbeginblock.write(code_file, instruction_number);
 
       if (linear_decomposition)
-        compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
-
-      // The equations
-      for (i = 0; i < block_size; i++)
-        {
-          //The Temporary terms
-          temporary_terms_t tt2;
-          if (v_temporary_terms[block][i].size() && !linear_decomposition)
-            {
-              for (auto it : v_temporary_terms[block][i])
-                {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
-                    it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
-
-                  FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(map_idx.find(it->idx)->second));
-                  fnumexpr.write(code_file, instruction_number);
-                  it->compile(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
-                  FSTPT_ fstpt(static_cast<int>(map_idx.find(it->idx)->second));
-                  fstpt.write(code_file, instruction_number);
-                  // Insert current node into tt2
-                  tt2.insert(it);
-#ifdef DEBUGC
-                  cout << "FSTPT " << v << endl;
-                  instruction_number++;
-                  code_file.write(&FOK, sizeof(FOK));
-                  code_file.write(reinterpret_cast<char *>(&k), sizeof(k));
-                  ki++;
-#endif
+        compileTemporaryTerms(code_file, instruction_number, true, false);
 
-                }
-            }
+      //The Temporary terms
+      deriv_node_temp_terms_t tef_terms;
+      if (!linear_decomposition)
+        for (auto it : blocks_temporary_terms[block])
+          {
+            if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+              it->compileExternalFunctionOutput(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false, tef_terms);
+
+            FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+            fnumexpr.write(code_file, instruction_number);
+            it->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false, tef_terms);
+            FSTPT_ fstpt(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+            fstpt.write(code_file, instruction_number);
 #ifdef DEBUGC
-          for (const auto &it : v_temporary_terms[block][i])
-            {
-              auto ii = map_idx.find(it->idx);
-              cout << "map_idx[" << it->idx <<"]=" << ii->second << endl;
-            }
+            cout << "FSTPT " << v << endl;
+            instruction_number++;
+            code_file.write(&FOK, sizeof(FOK));
+            code_file.write(reinterpret_cast<char *>(&k), sizeof(k));
+            ki++;
 #endif
+          }
 
+      // The equations
+      for (i = 0; i < block_size; i++)
+        {
           int variable_ID, equation_ID;
           EquationType equ_type;
 
@@ -1144,16 +1039,16 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
                   eq_node = getBlockEquationExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, true, false);
+                  rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
+                  lhs->compile(code_file, instruction_number, true, blocks_temporary_terms_idxs, true, false);
                 }
               else if (equ_type == EquationType::evaluate_s)
                 {
                   eq_node = getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, true, false);
+                  rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
+                  lhs->compile(code_file, instruction_number, true, blocks_temporary_terms_idxs, true, false);
                 }
               break;
             case BlockSimulationType::solveBackwardComplete:
@@ -1174,8 +1069,8 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
               eq_node = getBlockEquationExpr(block, i);
               lhs = eq_node->arg1;
               rhs = eq_node->arg2;
-              lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
-              rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+              lhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
+              rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
 
               FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
               fbinary.write(code_file, instruction_number);
@@ -1203,7 +1098,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
                 FNUMEXPR_ fnumexpr(FirstEndoDerivative, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0);
                 fnumexpr.write(code_file, instruction_number);
               }
-              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
+              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0);
               {
                 FSTPG_ fstpg(0);
                 fstpg.write(code_file, instruction_number);
@@ -1242,7 +1137,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
                       Uf[eqr].Ufl->lag = lag;
                       FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, lag);
                       fnumexpr.write(code_file, instruction_number);
-                      compileChainRuleDerivative(code_file, instruction_number, block, eq, var, lag, map_idx);
+                      compileChainRuleDerivative(code_file, instruction_number, block, eq, var, lag);
                       FSTPU_ fstpu(count_u);
                       fstpu.write(code_file, instruction_number);
                       count_u++;
@@ -1311,7 +1206,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
           int varr = getBlockVariableID(block, var);
           FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, lag);
           fnumexpr.write(code_file, instruction_number);
-          compileDerivative(code_file, instruction_number, eqr, varr, lag, map_idx);
+          compileDerivative(code_file, instruction_number, eqr, varr, lag);
           FSTPG3_ fstpg3(eq, var, lag, blocks_jacob_cols_endo[block].at({ var, lag }));
           fstpg3.write(code_file, instruction_number);
         }
@@ -1322,7 +1217,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
           int varr = 0; // Dummy value, actually unused by the bytecode MEX
           FNUMEXPR_ fnumexpr(FirstExoDerivative, eqr, varr, lag);
           fnumexpr.write(code_file, instruction_number);
-          d->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+          d->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
           FSTPG3_ fstpg3(eq, var, lag, blocks_jacob_cols_exo[block].at({ var, lag }));
           fstpg3.write(code_file, instruction_number);
         }
@@ -1333,7 +1228,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
           int varr = 0; // Dummy value, actually unused by the bytecode MEX
           FNUMEXPR_ fnumexpr(FirstExodetDerivative, eqr, varr, lag);
           fnumexpr.write(code_file, instruction_number);
-          d->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+          d->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
           FSTPG3_ fstpg3(eq, var, lag, blocks_jacob_cols_exo_det[block].at({ var, lag }));
           fstpg3.write(code_file, instruction_number);
         }
@@ -1344,7 +1239,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
           int varr = 0; // Dummy value, actually unused by the bytecode MEX
           FNUMEXPR_ fnumexpr(FirstOtherEndoDerivative, eqr, varr, lag);
           fnumexpr.write(code_file, instruction_number);
-          d->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
+          d->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, true, false);
           FSTPG3_ fstpg3(eq, var, lag, blocks_jacob_cols_other_endo[block].at({ var, lag }));
           fstpg3.write(code_file, instruction_number);
         }
@@ -1648,31 +1543,29 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
   mDynamicModelFile << "function [varargout] = dynamic(options_, M_, oo_, varargin)" << endl
                     << "  g2=[];g3=[];" << endl;
   //Temporary variables declaration
-  bool OK = true;
-  ostringstream tmp_output;
-  for (auto temporary_term : temporary_terms)
+  if (blocks_temporary_terms_idxs.size() > 0)
     {
-      if (OK)
-        OK = false;
-      else
-        tmp_output << " ";
-      // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-      temporary_term->writeOutput(tmp_output, ExprNodeOutputType::matlabStaticModelSparse, temporary_terms, {});
-    }
-  if (tmp_output.str().length() > 0)
-    mDynamicModelFile << "  global " << tmp_output.str() << ";" << endl;
+      temporary_terms_t tt_all;
+      for (auto [tt, idx] : blocks_temporary_terms_idxs)
+        tt_all.insert(tt);
 
-  mDynamicModelFile << "  T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);" << endl;
-  tmp_output.str("");
-  for (auto temporary_term : temporary_terms)
-    {
-      tmp_output << "  ";
-      // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-      temporary_term->writeOutput(tmp_output, ExprNodeOutputType::matlabStaticModelSparse, temporary_terms, {});
-      tmp_output << "=T_init;" << endl;
+      mDynamicModelFile << "  global";
+      for (auto tt : tt_all)
+        {
+          mDynamicModelFile << " ";
+          // In the following, "Static" is used to avoid getting the "(it_)" subscripting
+          tt->writeOutput(mDynamicModelFile, ExprNodeOutputType::matlabStaticModelSparse, tt_all, {});
+        }
+      mDynamicModelFile << ";" << endl
+                        << "  T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);" << endl;
+      for (auto tt : tt_all)
+        {
+          mDynamicModelFile << "  ";
+          // In the following, "Static" is used to avoid getting the "(it_)" subscripting
+          tt->writeOutput(mDynamicModelFile, ExprNodeOutputType::matlabStaticModelSparse, tt_all, {});
+          mDynamicModelFile << "=T_init;" << endl;
+        }
     }
-  if (tmp_output.str().length() > 0)
-    mDynamicModelFile << tmp_output.str();
 
   mDynamicModelFile << "  y_kmin=M_.maximum_lag;" << endl
                     << "  y_kmax=M_.maximum_lead;" << endl
@@ -4525,9 +4418,8 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
 
       computeBlockDynJacobianCols();
 
-      global_temporary_terms = true;
       if (!no_tmp_terms)
-        computeTemporaryTermsOrdered();
+        computeBlockTemporaryTerms();
     }
   else if (block)
     {
@@ -4555,15 +4447,12 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
 
       computeBlockDynJacobianCols();
 
-      global_temporary_terms = true;
       if (!no_tmp_terms)
-        computeTemporaryTermsOrdered();
+        computeBlockTemporaryTerms();
     }
   else
     {
       computeTemporaryTerms(!use_dll, no_tmp_terms);
-      if (bytecode && !no_tmp_terms)
-        computeTemporaryTermsMapping();
 
       /* Must be called after computeTemporaryTerms(), because it depends on
          temporary_terms_mlv to be filled */
@@ -4832,12 +4721,12 @@ void
 DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const
 {
   if (block && bytecode)
-    writeModelEquationsCode_Block(basename, map_idx, linear_decomposition);
+    writeModelEquationsCode_Block(basename, linear_decomposition);
   else if (!block && bytecode)
     {
       if (linear_decomposition)
-        writeModelEquationsCode_Block(basename, map_idx, linear_decomposition);
-      writeModelEquationsCode(basename, map_idx);
+        writeModelEquationsCode_Block(basename, linear_decomposition);
+      writeModelEquationsCode(basename);
     }
   else if (block && !bytecode)
     writeSparseDynamicMFile(basename);
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index c61c2e9e21389b6ca9706a2e04e4359f55d77dc8..75e815436e01d4c31956aa5228ffa3df21442984 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -136,9 +136,9 @@ private:
   //! Writes the Block reordred structure of the model in M output
   void writeModelEquationsOrdered_M(const string &basename) const;
   //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
-  void writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx, bool linear_decomposition) const;
+  void writeModelEquationsCode_Block(const string &basename, bool linear_decomposition) const;
   //! Writes the code of the model in virtual machine bytecode
-  void writeModelEquationsCode(const string &basename, const map_idx_t &map_idx) const;
+  void writeModelEquationsCode(const string &basename) const;
 
   void writeSetAuxiliaryVariables(const string &basename, bool julia) const;
   void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const;
@@ -161,15 +161,14 @@ private:
 
   string reform(const string &name) const;
 
-  //! sorts the temporary terms in the blocks order
-  void computeTemporaryTermsOrdered();
+  void additionalBlockTemporaryTerms(int blk,
+                                     vector<temporary_terms_t> &blocks_temporary_terms,
+                                     map<expr_t, pair<int, int>> &reference_count) const override;
 
-  //! creates a mapping from the index of temporary terms to a natural index
-  void computeTemporaryTermsMapping();
   //! Write derivative code of an equation w.r. to a variable
-  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const;
+  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag) const;
   //! Write chain rule derivative code of an equation w.r. to a variable
-  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const map_idx_t &map_idx) const;
+  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag) const;
 
   //! Get the type corresponding to a derivation ID
   SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
@@ -199,9 +198,6 @@ private:
   */
   void substituteLeadLagInternal(AuxVarType type, bool deterministic_model, const vector<string> &subset);
 
-  //! Indicate if the temporary terms are computed for the overall model (true) or not (false). Default value true
-  bool global_temporary_terms{true};
-
   //! Help computeXrefs to compute the reverse references (i.e. param->eqs, endo->eqs, etc)
   void computeRevXref(map<pair<int, int>, set<int>> &xrefset, const set<pair<int, int>> &eiref, int eqn);
 
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 7de1de6854f8a2cc9d8b5a601fcab3403e96d129..74d741fe2be484c65fde0623d8eab1a8f0c7715b 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -75,7 +75,7 @@ ExprNode::cost(int cost, bool is_matlab) const
 }
 
 int
-ExprNode::cost(const temporary_terms_t &temp_terms_map, bool is_matlab) const
+ExprNode::cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const
 {
   // For a terminal node, the cost is null
   return 0;
@@ -161,12 +161,8 @@ ExprNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
 }
 
 void
-ExprNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                temporary_terms_t &temporary_terms,
-                                map<expr_t, pair<int, int>> &first_occurence,
-                                int Curr_block,
-                                vector<vector<temporary_terms_t>> &v_temporary_terms,
-                                int equation) const
+ExprNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                     map<expr_t, pair<int, int>> &reference_count) const
 {
   // Nothing to do for a terminal node
 }
@@ -191,10 +187,10 @@ ExprNode::writeOutput(ostream &output, ExprNodeOutputType output_type, const tem
 
 void
 ExprNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                  bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                  const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+                  bool lhs_rhs,
+                  const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic) const
 {
-  compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, {});
+  compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, steady_dynamic, {});
 }
 
 void
@@ -217,8 +213,8 @@ ExprNode::writeJsonExternalFunctionOutput(vector<string> &efout,
 
 void
 ExprNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                        const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                        bool lhs_rhs,
+                                        const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                         deriv_node_temp_terms_t &tef_terms) const
 {
   // Nothing to do
@@ -423,13 +419,6 @@ NumConstNode::computeDerivative(int deriv_id)
   return datatree.Zero;
 }
 
-void
-NumConstNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
-{
-  if (temporary_terms.find(const_cast<NumConstNode *>(this)) != temporary_terms.end())
-    temporary_terms_inuse.insert(idx);
-}
-
 void
 NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                           const temporary_terms_t &temporary_terms,
@@ -470,8 +459,8 @@ NumConstNode::eval(const eval_context_t &eval_context) const noexcept(false)
 
 void
 NumConstNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                      bool lhs_rhs,
+                      const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                       const deriv_node_temp_terms_t &tef_terms) const
 {
   FLDC_ fldc(datatree.num_constants.getDouble(id));
@@ -873,15 +862,6 @@ VariableNode::computeDerivative(int deriv_id)
   exit(EXIT_FAILURE);
 }
 
-void
-VariableNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
-{
-  if (temporary_terms.find(const_cast<VariableNode *>(this)) != temporary_terms.end())
-    temporary_terms_inuse.insert(idx);
-  if (get_type() == SymbolType::modelLocalVariable)
-    datatree.getLocalVariable(symb_id)->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
-}
-
 bool
 VariableNode::containsExternalFunction() const
 {
@@ -1252,13 +1232,13 @@ VariableNode::eval(const eval_context_t &eval_context) const noexcept(false)
 
 void
 VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                      bool lhs_rhs,
+                      const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                       const deriv_node_temp_terms_t &tef_terms) const
 {
   auto type = get_type();
   if (type == SymbolType::modelLocalVariable || type == SymbolType::modFileLocalVariable)
-    datatree.getLocalVariable(symb_id)->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
+    datatree.getLocalVariable(symb_id)->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
   else
     {
       int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
@@ -1325,18 +1305,6 @@ VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number,
     }
 }
 
-void
-VariableNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                    temporary_terms_t &temporary_terms,
-                                    map<expr_t, pair<int, int>> &first_occurence,
-                                    int Curr_block,
-                                    vector<vector<temporary_terms_t>> &v_temporary_terms,
-                                    int equation) const
-{
-  if (get_type() == SymbolType::modelLocalVariable)
-    datatree.getLocalVariable(symb_id)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
-}
-
 void
 VariableNode::collectVARLHSVariable(set<expr_t> &result) const
 {
@@ -2224,13 +2192,14 @@ UnaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map,
 }
 
 int
-UnaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) const
+UnaryOpNode::cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const
 {
   // For a temporary term, the cost is null
-  if (temporary_terms.find(const_cast<UnaryOpNode *>(this)) != temporary_terms.end())
-    return 0;
+  for (const auto &it : blocks_temporary_terms)
+    if (it.find(const_cast<UnaryOpNode *>(this)) != it.end())
+      return 0;
 
-  return cost(arg->cost(temporary_terms, is_matlab), is_matlab);
+  return cost(arg->cost(blocks_temporary_terms, is_matlab), is_matlab);
 }
 
 int
@@ -2362,41 +2331,25 @@ UnaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
 }
 
 void
-UnaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                   temporary_terms_t &temporary_terms,
-                                   map<expr_t, pair<int, int>> &first_occurence,
-                                   int Curr_block,
-                                   vector< vector<temporary_terms_t>> &v_temporary_terms,
-                                   int equation) const
+UnaryOpNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                        map<expr_t, pair<int, int>> &reference_count) const
 {
   expr_t this2 = const_cast<UnaryOpNode *>(this);
   if (auto it = reference_count.find(this2);
       it == reference_count.end())
     {
-      reference_count[this2] = 1;
-      first_occurence[this2] = { Curr_block, equation };
-      arg->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
+      reference_count[this2] = { 1, blk };
+      arg->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
     }
   else
     {
-      reference_count[this2]++;
-      if (reference_count[this2] * cost(temporary_terms, false) > min_cost_c)
-        {
-          temporary_terms.insert(this2);
-          v_temporary_terms[first_occurence[this2].first][first_occurence[this2].second].insert(this2);
-        }
+      auto &[nref, first_blk] = it->second;
+      nref++;
+      if (nref * cost(blocks_temporary_terms, false) > min_cost_c)
+        blocks_temporary_terms[first_blk].insert(this2);
     }
 }
 
-void
-UnaryOpNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
-{
-  if (temporary_terms.find(const_cast<UnaryOpNode *>(this)) != temporary_terms.end())
-    temporary_terms_inuse.insert(idx);
-  else
-    arg->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
-}
-
 bool
 UnaryOpNode::containsExternalFunction() const
 {
@@ -2912,11 +2865,11 @@ UnaryOpNode::writeJsonExternalFunctionOutput(vector<string> &efout,
 
 void
 UnaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                           const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                           bool lhs_rhs,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                            deriv_node_temp_terms_t &tef_terms) const
 {
-  arg->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+  arg->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs,
                                      dynamic, steady_dynamic, tef_terms);
 }
 
@@ -2998,31 +2951,30 @@ UnaryOpNode::eval(const eval_context_t &eval_context) const noexcept(false)
 
 void
 UnaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                     const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                     bool lhs_rhs,
+                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                      const deriv_node_temp_terms_t &tef_terms) const
 {
-  if (temporary_terms.find(const_cast<UnaryOpNode *>(this)) != temporary_terms.end())
+  if (auto it = temporary_terms_idxs.find(const_cast<UnaryOpNode *>(this));
+      it != temporary_terms_idxs.end())
     {
       if (dynamic)
         {
-          auto ii = map_idx.find(idx);
-          FLDT_ fldt(ii->second);
+          FLDT_ fldt(it->second);
           fldt.write(CompileCode, instruction_number);
         }
       else
         {
-          auto ii = map_idx.find(idx);
-          FLDST_ fldst(ii->second);
+          FLDST_ fldst(it->second);
           fldst.write(CompileCode, instruction_number);
         }
       return;
     }
   if (op_code == UnaryOpcode::steadyState)
-    arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, true, tef_terms);
+    arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, true, tef_terms);
   else
     {
-      arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
+      arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
       FUNARY_ funary{static_cast<uint8_t>(op_code)};
       funary.write(CompileCode, instruction_number);
     }
@@ -4089,13 +4041,14 @@ BinaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map,
 }
 
 int
-BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) const
+BinaryOpNode::cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const
 {
   // For a temporary term, the cost is null
-  if (temporary_terms.find(const_cast<BinaryOpNode *>(this)) != temporary_terms.end())
-    return 0;
+  for (const auto &it : blocks_temporary_terms)
+    if (it.find(const_cast<BinaryOpNode *>(this)) != it.end())
+      return 0;
 
-  int arg_cost = arg1->cost(temporary_terms, is_matlab) + arg2->cost(temporary_terms, is_matlab);
+  int arg_cost = arg1->cost(blocks_temporary_terms, is_matlab) + arg2->cost(blocks_temporary_terms, is_matlab);
 
   return cost(arg_cost, is_matlab);
 }
@@ -4189,31 +4142,24 @@ BinaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
 }
 
 void
-BinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                    temporary_terms_t &temporary_terms,
-                                    map<expr_t, pair<int, int>> &first_occurence,
-                                    int Curr_block,
-                                    vector<vector<temporary_terms_t>> &v_temporary_terms,
-                                    int equation) const
+BinaryOpNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                         map<expr_t, pair<int, int>> &reference_count) const
 {
   expr_t this2 = const_cast<BinaryOpNode *>(this);
   if (auto it = reference_count.find(this2);
       it == reference_count.end())
     {
-      reference_count[this2] = 1;
-      first_occurence[this2] = { Curr_block, equation };
-      arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
-      arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
+      reference_count[this2] = { 1, blk };
+      arg1->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+      arg2->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
     }
   else
     {
-      reference_count[this2]++;
-      if (reference_count[this2] * cost(temporary_terms, false) > min_cost_c
+      auto &[nref, first_blk] = it->second;
+      nref++;
+      if (nref * cost(blocks_temporary_terms, false) > min_cost_c
           && op_code != BinaryOpcode::equal)
-        {
-          temporary_terms.insert(this2);
-          v_temporary_terms[first_occurence[this2].first][first_occurence[this2].second].insert(this2);
-        }
+        blocks_temporary_terms[first_blk].insert(this2);
     }
 }
 
@@ -4283,23 +4229,22 @@ BinaryOpNode::eval(const eval_context_t &eval_context) const noexcept(false)
 
 void
 BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                      bool lhs_rhs,
+                      const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                       const deriv_node_temp_terms_t &tef_terms) const
 {
   // If current node is a temporary term
-  if (temporary_terms.find(const_cast<BinaryOpNode *>(this)) != temporary_terms.end())
+  if (auto it = temporary_terms_idxs.find(const_cast<BinaryOpNode *>(this));
+      it != temporary_terms_idxs.end())
     {
       if (dynamic)
         {
-          auto ii = map_idx.find(idx);
-          FLDT_ fldt(ii->second);
+          FLDT_ fldt(it->second);
           fldt.write(CompileCode, instruction_number);
         }
       else
         {
-          auto ii = map_idx.find(idx);
-          FLDST_ fldst(ii->second);
+          FLDST_ fldst(it->second);
           fldst.write(CompileCode, instruction_number);
         }
       return;
@@ -4309,24 +4254,12 @@ BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
       FLDC_ fldc(powerDerivOrder);
       fldc.write(CompileCode, instruction_number);
     }
-  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
-  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
+  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
   FBINARY_ fbinary{static_cast<int>(op_code)};
   fbinary.write(CompileCode, instruction_number);
 }
 
-void
-BinaryOpNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
-{
-  if (temporary_terms.find(const_cast<BinaryOpNode *>(this)) != temporary_terms.end())
-    temporary_terms_inuse.insert(idx);
-  else
-    {
-      arg1->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
-      arg2->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
-    }
-}
-
 bool
 BinaryOpNode::containsExternalFunction() const
 {
@@ -4718,13 +4651,13 @@ BinaryOpNode::writeJsonExternalFunctionOutput(vector<string> &efout,
 
 void
 BinaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                            bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                            const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                            bool lhs_rhs,
+                                            const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                             deriv_node_temp_terms_t &tef_terms) const
 {
-  arg1->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+  arg1->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs,
                                       dynamic, steady_dynamic, tef_terms);
-  arg2->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+  arg2->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs,
                                       dynamic, steady_dynamic, tef_terms);
 }
 
@@ -5905,15 +5838,16 @@ TrinaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map
 }
 
 int
-TrinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) const
+TrinaryOpNode::cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const
 {
   // For a temporary term, the cost is null
-  if (temporary_terms.find(const_cast<TrinaryOpNode *>(this)) != temporary_terms.end())
-    return 0;
+  for (const auto &it : blocks_temporary_terms)
+    if (it.find(const_cast<TrinaryOpNode *>(this)) != it.end())
+      return 0;
 
-  int arg_cost = arg1->cost(temporary_terms, is_matlab)
-    + arg2->cost(temporary_terms, is_matlab)
-    + arg3->cost(temporary_terms, is_matlab);
+  int arg_cost = arg1->cost(blocks_temporary_terms, is_matlab)
+    + arg2->cost(blocks_temporary_terms, is_matlab)
+    + arg3->cost(blocks_temporary_terms, is_matlab);
 
   return cost(arg_cost, is_matlab);
 }
@@ -5969,31 +5903,24 @@ TrinaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
 }
 
 void
-TrinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                     temporary_terms_t &temporary_terms,
-                                     map<expr_t, pair<int, int>> &first_occurence,
-                                     int Curr_block,
-                                     vector<vector<temporary_terms_t>> &v_temporary_terms,
-                                     int equation) const
+TrinaryOpNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                          map<expr_t, pair<int, int>> &reference_count) const
 {
   expr_t this2 = const_cast<TrinaryOpNode *>(this);
   if (auto it = reference_count.find(this2);
       it == reference_count.end())
     {
-      reference_count[this2] = 1;
-      first_occurence[this2] = { Curr_block, equation };
-      arg1->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
-      arg2->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
-      arg3->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, Curr_block, v_temporary_terms, equation);
+      reference_count[this2] = { 1, blk };
+      arg1->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+      arg2->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+      arg3->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
     }
   else
     {
-      reference_count[this2]++;
-      if (reference_count[this2] * cost(temporary_terms, false) > min_cost_c)
-        {
-          temporary_terms.insert(this2);
-          v_temporary_terms[first_occurence[this2].first][first_occurence[this2].second].insert(this2);
-        }
+      auto &[nref, first_blk] = it->second;
+      nref++;
+      if (nref * cost(blocks_temporary_terms, false) > min_cost_c)
+        blocks_temporary_terms[first_blk].insert(this2);
     }
 }
 
@@ -6023,47 +5950,33 @@ TrinaryOpNode::eval(const eval_context_t &eval_context) const noexcept(false)
 
 void
 TrinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                       bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                       const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                       bool lhs_rhs,
+                       const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                        const deriv_node_temp_terms_t &tef_terms) const
 {
   // If current node is a temporary term
-  if (temporary_terms.find(const_cast<TrinaryOpNode *>(this)) != temporary_terms.end())
+  if (auto it = temporary_terms_idxs.find(const_cast<TrinaryOpNode *>(this));
+      it != temporary_terms_idxs.end())
     {
       if (dynamic)
         {
-          auto ii = map_idx.find(idx);
-          FLDT_ fldt(ii->second);
+          FLDT_ fldt(it->second);
           fldt.write(CompileCode, instruction_number);
         }
       else
         {
-          auto ii = map_idx.find(idx);
-          FLDST_ fldst(ii->second);
+          FLDST_ fldst(it->second);
           fldst.write(CompileCode, instruction_number);
         }
       return;
     }
-  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
-  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
-  arg3->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
+  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg3->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
   FTRINARY_ ftrinary{static_cast<int>(op_code)};
   ftrinary.write(CompileCode, instruction_number);
 }
 
-void
-TrinaryOpNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
-{
-  if (temporary_terms.find(const_cast<TrinaryOpNode *>(this)) != temporary_terms.end())
-    temporary_terms_inuse.insert(idx);
-  else
-    {
-      arg1->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
-      arg2->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
-      arg3->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
-    }
-}
-
 bool
 TrinaryOpNode::containsExternalFunction() const
 {
@@ -6212,15 +6125,15 @@ TrinaryOpNode::writeJsonExternalFunctionOutput(vector<string> &efout,
 
 void
 TrinaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                             bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                             const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                             bool lhs_rhs,
+                                             const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                              deriv_node_temp_terms_t &tef_terms) const
 {
-  arg1->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+  arg1->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs,
                                       dynamic, steady_dynamic, tef_terms);
-  arg2->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+  arg2->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs,
                                       dynamic, steady_dynamic, tef_terms);
-  arg3->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+  arg3->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs,
                                       dynamic, steady_dynamic, tef_terms);
 }
 
@@ -6720,12 +6633,12 @@ AbstractExternalFunctionNode::getChainRuleDerivative(int deriv_id, const map<int
 
 unsigned int
 AbstractExternalFunctionNode::compileExternalFunctionArguments(ostream &CompileCode, unsigned int &instruction_number,
-                                                               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                               const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                               bool lhs_rhs,
+                                                               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                                                const deriv_node_temp_terms_t &tef_terms) const
 {
   for (auto argument : arguments)
-    argument->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+    argument->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms_idxs,
                       dynamic, steady_dynamic, tef_terms);
   return (arguments.size());
 }
@@ -6744,16 +6657,6 @@ AbstractExternalFunctionNode::collectDynamicVariables(SymbolType type_arg, set<p
     argument->collectDynamicVariables(type_arg, result);
 }
 
-void
-AbstractExternalFunctionNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
-{
-  if (temporary_terms.find(const_cast<AbstractExternalFunctionNode *>(this)) != temporary_terms.end())
-    temporary_terms_inuse.insert(idx);
-  else
-    for (auto argument : arguments)
-      argument->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
-}
-
 double
 AbstractExternalFunctionNode::eval(const eval_context_t &eval_context) const noexcept(false)
 {
@@ -7068,6 +6971,23 @@ AbstractExternalFunctionNode::computeTemporaryTerms(const pair<int, int> &derivO
   temp_terms_map[derivOrder].insert(const_cast<AbstractExternalFunctionNode *>(this));
 }
 
+void
+AbstractExternalFunctionNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                                         map<expr_t, pair<int, int>> &reference_count) const
+{
+  // See comments in computeTemporaryTerms() for the logic
+  expr_t this2 = const_cast<AbstractExternalFunctionNode *>(this);
+  for (auto &tt : blocks_temporary_terms)
+    if (auto it = find_if(tt.cbegin(), tt.cend(), sameTefTermPredicate());
+        it != tt.cend())
+      {
+        tt.insert(this2);
+        return;
+      }
+
+  blocks_temporary_terms[blk].insert(this2);
+}
+
 bool
 AbstractExternalFunctionNode::isNumConstNodeEqualTo(double value) const
 {
@@ -7300,38 +7220,23 @@ ExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
   return theDeriv;
 }
 
-void
-ExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                            temporary_terms_t &temporary_terms,
-                                            map<expr_t, pair<int, int>> &first_occurence,
-                                            int Curr_block,
-                                            vector< vector<temporary_terms_t>> &v_temporary_terms,
-                                            int equation) const
-{
-  expr_t this2 = const_cast<ExternalFunctionNode *>(this);
-  temporary_terms.insert(this2);
-  first_occurence[this2] = { Curr_block, equation };
-  v_temporary_terms[Curr_block][equation].insert(this2);
-}
-
 void
 ExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                              bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                              const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                              bool lhs_rhs,
+                              const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                               const deriv_node_temp_terms_t &tef_terms) const
 {
-  if (temporary_terms.find(const_cast<ExternalFunctionNode *>(this)) != temporary_terms.end())
+  if (auto it = temporary_terms_idxs.find(const_cast<ExternalFunctionNode *>(this));
+      it != temporary_terms_idxs.end())
     {
       if (dynamic)
         {
-          auto ii = map_idx.find(idx);
-          FLDT_ fldt(ii->second);
+          FLDT_ fldt(it->second);
           fldt.write(CompileCode, instruction_number);
         }
       else
         {
-          auto ii = map_idx.find(idx);
-          FLDST_ fldst(ii->second);
+          FLDST_ fldst(it->second);
           fldst.write(CompileCode, instruction_number);
         }
       return;
@@ -7351,16 +7256,16 @@ ExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_nu
 
 void
 ExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    bool lhs_rhs,
+                                                    const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                                     deriv_node_temp_terms_t &tef_terms) const
 {
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
   assert(first_deriv_symb_id != ExternalFunctionsTable::IDSetButNoNameProvided);
 
   for (auto argument : arguments)
-    argument->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                            map_idx, dynamic, steady_dynamic, tef_terms);
+    argument->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs,
+                                            temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 
   if (!alreadyWrittenAsTefTerm(symb_id, tef_terms))
     {
@@ -7377,8 +7282,8 @@ ExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsign
         nb_output_arguments = 2;
       else
         nb_output_arguments = 1;
-      unsigned int nb_input_arguments = compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                                                         map_idx, dynamic, steady_dynamic, tef_terms);
+      unsigned int nb_input_arguments = compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs,
+                                                                         temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 
       FCALL_ fcall(nb_output_arguments, nb_input_arguments, datatree.symbol_table.getName(symb_id), indx);
       switch (nb_output_arguments)
@@ -7616,20 +7521,6 @@ FirstDerivExternalFunctionNode::FirstDerivExternalFunctionNode(DataTree &datatre
 {
 }
 
-void
-FirstDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                                      temporary_terms_t &temporary_terms,
-                                                      map<expr_t, pair<int, int>> &first_occurence,
-                                                      int Curr_block,
-                                                      vector< vector<temporary_terms_t>> &v_temporary_terms,
-                                                      int equation) const
-{
-  expr_t this2 = const_cast<FirstDerivExternalFunctionNode *>(this);
-  temporary_terms.insert(this2);
-  first_occurence[this2] = { Curr_block, equation };
-  v_temporary_terms[Curr_block][equation].insert(this2);
-}
-
 expr_t
 FirstDerivExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 {
@@ -7721,22 +7612,21 @@ FirstDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
 
 void
 FirstDerivExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                        const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                        bool lhs_rhs,
+                                        const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                         const deriv_node_temp_terms_t &tef_terms) const
 {
-  if (temporary_terms.find(const_cast<FirstDerivExternalFunctionNode *>(this)) != temporary_terms.end())
+  if (auto it = temporary_terms_idxs.find(const_cast<FirstDerivExternalFunctionNode *>(this));
+      it != temporary_terms_idxs.end())
     {
       if (dynamic)
         {
-          auto ii = map_idx.find(idx);
-          FLDT_ fldt(ii->second);
+          FLDT_ fldt(it->second);
           fldt.write(CompileCode, instruction_number);
         }
       else
         {
-          auto ii = map_idx.find(idx);
-          FLDST_ fldst(ii->second);
+          FLDST_ fldst(it->second);
           fldst.write(CompileCode, instruction_number);
         }
       return;
@@ -7903,8 +7793,8 @@ FirstDerivExternalFunctionNode::writeJsonExternalFunctionOutput(vector<string> &
 
 void
 FirstDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                                              bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                              const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                              bool lhs_rhs,
+                                                              const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                                               deriv_node_temp_terms_t &tef_terms) const
 {
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
@@ -7913,8 +7803,8 @@ FirstDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCo
   if (first_deriv_symb_id == symb_id || alreadyWrittenAsTefTerm(first_deriv_symb_id, tef_terms))
     return;
 
-  unsigned int nb_add_input_arguments = compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                                                         map_idx, dynamic, steady_dynamic, tef_terms);
+  unsigned int nb_add_input_arguments = compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs,
+                                                                         temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
   if (first_deriv_symb_id == ExternalFunctionsTable::IDNotSet)
     {
       unsigned int nb_input_arguments = 0;
@@ -8008,20 +7898,6 @@ SecondDerivExternalFunctionNode::SecondDerivExternalFunctionNode(DataTree &datat
 {
 }
 
-void
-SecondDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                                       temporary_terms_t &temporary_terms,
-                                                       map<expr_t, pair<int, int>> &first_occurence,
-                                                       int Curr_block,
-                                                       vector< vector<temporary_terms_t>> &v_temporary_terms,
-                                                       int equation) const
-{
-  expr_t this2 = const_cast<SecondDerivExternalFunctionNode *>(this);
-  temporary_terms.insert(this2);
-  first_occurence[this2] = { Curr_block, equation };
-  v_temporary_terms[Curr_block][equation].insert(this2);
-}
-
 expr_t
 SecondDerivExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 {
@@ -8303,8 +8179,8 @@ SecondDerivExternalFunctionNode::computeXrefs(EquationInfo &ei) const
 
 void
 SecondDerivExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                                         bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                         const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                         bool lhs_rhs,
+                                         const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                          const deriv_node_temp_terms_t &tef_terms) const
 {
   cerr << "SecondDerivExternalFunctionNode::compile: not implemented." << endl;
@@ -8313,8 +8189,8 @@ SecondDerivExternalFunctionNode::compile(ostream &CompileCode, unsigned int &ins
 
 void
 SecondDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                                               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                               const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                               bool lhs_rhs,
+                                                               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                                                deriv_node_temp_terms_t &tef_terms) const
 {
   cerr << "SecondDerivExternalFunctionNode::compileExternalFunctionOutput: not implemented." << endl;
@@ -8356,14 +8232,10 @@ VarExpectationNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
 }
 
 void
-VarExpectationNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                          temporary_terms_t &temporary_terms,
-                                          map<expr_t, pair<int, int>> &first_occurence,
-                                          int Curr_block,
-                                          vector< vector<temporary_terms_t>> &v_temporary_terms,
-                                          int equation) const
+VarExpectationNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                               map<expr_t, pair<int, int>> &reference_count) const
 {
-  cerr << "VarExpectationNode::computeTemporaryTerms not implemented." << endl;
+  cerr << "VarExpectationNode::computeBlocksTemporaryTerms not implemented." << endl;
   exit(EXIT_FAILURE);
 }
 
@@ -8551,17 +8423,10 @@ VarExpectationNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, i
 {
 }
 
-void
-VarExpectationNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
-{
-  cerr << "VarExpectationNode::collectTemporary_terms not implemented." << endl;
-  exit(EXIT_FAILURE);
-}
-
 void
 VarExpectationNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                            bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                            const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                            bool lhs_rhs,
+                            const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                             const deriv_node_temp_terms_t &tef_terms) const
 {
   cerr << "VarExpectationNode::compile not implemented." << endl;
@@ -8813,17 +8678,10 @@ PacExpectationNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
 }
 
 void
-PacExpectationNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                          temporary_terms_t &temporary_terms,
-                                          map<expr_t, pair<int, int>> &first_occurence,
-                                          int Curr_block,
-                                          vector< vector<temporary_terms_t>> &v_temporary_terms,
-                                          int equation) const
+PacExpectationNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                               map<expr_t, pair<int, int>> &reference_count) const
 {
-  expr_t this2 = const_cast<PacExpectationNode *>(this);
-  temporary_terms.insert(this2);
-  first_occurence[this2] = { Curr_block, equation };
-  v_temporary_terms[Curr_block][equation].insert(this2);
+  blocks_temporary_terms[blk].insert(const_cast<PacExpectationNode *>(this));
 }
 
 expr_t
@@ -8981,17 +8839,10 @@ PacExpectationNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, i
 {
 }
 
-void
-PacExpectationNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
-{
-  if (temporary_terms.find(const_cast<PacExpectationNode *>(this)) != temporary_terms.end())
-    temporary_terms_inuse.insert(idx);
-}
-
 void
 PacExpectationNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                            bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                            const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                            bool lhs_rhs,
+                            const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                             const deriv_node_temp_terms_t &tef_terms) const
 {
   cerr << "PacExpectationNode::compile not implemented." << endl;
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 9f8bcd15bc2cf71ec35d3591458299257b616043..e199a3e3371ffed0b64339bc7fdc99b20ec15df9 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2019 Dynare Team
+ * Copyright © 2007-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -49,11 +49,6 @@ using temporary_terms_t = set<expr_t, ExprNodeLess>;
 /*! Keeps track of array indices of temporary_terms for writing */
 using temporary_terms_idxs_t = map<expr_t, int>;
 
-//! set of temporary terms used in a block
-using temporary_terms_inuse_t = set<int>;
-
-using map_idx_t = map<int, int>;
-
 //! Type for evaluation contexts
 /*! The key is a symbol id. Lags are assumed to be null */
 using eval_context_t = map<int, double>;
@@ -230,7 +225,7 @@ protected:
   //! Cost of computing current node
   /*! Nodes included in temporary_terms are considered having a null cost */
   virtual int cost(int cost, bool is_matlab) const;
-  virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual int cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const;
   virtual int cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const;
 
   //! For creating equation cross references
@@ -306,6 +301,21 @@ public:
                                      map<expr_t, pair<int, pair<int, int>>> &reference_count,
                                      bool is_matlab) const;
 
+  //! Compute temporary terms in this expression for block decomposed model
+  /*!
+    \param[in] blk the block number
+    \param[out] blocks_temporary_terms the computed temporary terms, per block
+    \param[out] reference_count a temporary structure, used to count
+    references to each node (first integer is the
+    reference count, second integer is the number of the block in which the
+    expression first appears)
+
+    Same rules as computeTemporaryTerms() for computing cost.
+  */
+  virtual void computeBlockTemporaryTerms(int blk,
+                                          vector<temporary_terms_t> &blocks_temporary_terms,
+                                          map<expr_t, pair<int, int>> &reference_count) const;
+
   //! Writes output of node, using a Txxx notation for nodes in temporary_terms, and specifiying the set of already written external functions
   /*!
     \param[in] output the output stream
@@ -352,8 +362,8 @@ public:
                                                bool isdynamic = true) const;
 
   virtual void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                             bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                             const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                             bool lhs_rhs,
+                                             const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                              deriv_node_temp_terms_t &tef_terms) const;
 
   //! Computes the set of all variables of a given symbol type in the expression (with information on lags)
@@ -399,15 +409,6 @@ public:
   */
   virtual void collectExogenous(set<pair<int, int>> &result) const;
 
-  virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const = 0;
-
-  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                                     temporary_terms_t &temporary_terms,
-                                     map<expr_t, pair<int, int>> &first_occurence,
-                                     int Curr_block,
-                                     vector< vector<temporary_terms_t>> &v_temporary_terms,
-                                     int equation) const;
-
   class EvalException
   {
   };
@@ -417,8 +418,8 @@ public:
   };
 
   virtual double eval(const eval_context_t &eval_context) const noexcept(false) = 0;
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const = 0;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const = 0;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic) const;
   //! Creates a static version of this node
   /*!
     This method duplicates the current node by creating a similar node from which all leads/lags have been stripped,
@@ -749,9 +750,8 @@ public:
   bool containsExternalFunction() const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
-  void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const override;
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
@@ -826,15 +826,8 @@ public:
   bool containsExternalFunction() const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
-  void computeTemporaryTerms(map<expr_t, int > &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
-  void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const override;
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   SymbolType get_type() const;
@@ -906,7 +899,7 @@ public:
 private:
   expr_t computeDerivative(int deriv_id) override;
   int cost(int cost, bool is_matlab) const override;
-  int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const override;
+  int cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const override;
   int cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const override;
   //! Returns the derivative of this node if darg is the derivative of the argument
   expr_t composeDerivatives(expr_t darg, int deriv_id);
@@ -917,6 +910,8 @@ public:
                              map<pair<int, int>, temporary_terms_t> &temp_terms_map,
                              map<expr_t, pair<int, pair<int, int>>> &reference_count,
                              bool is_matlab) const override;
+  void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                  map<expr_t, pair<int, int>> &reference_count) const override;
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
@@ -930,21 +925,14 @@ public:
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
   void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                     bool lhs_rhs,
+                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                      deriv_node_temp_terms_t &tef_terms) const override;
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
-  void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const override;
   static double eval_opcode(UnaryOpcode op_code, double v) noexcept(false);
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
@@ -1015,7 +1003,7 @@ public:
 private:
   expr_t computeDerivative(int deriv_id) override;
   int cost(int cost, bool is_matlab) const override;
-  int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const override;
+  int cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const override;
   int cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const override;
   //! Returns the derivative of this node if darg1 and darg2 are the derivatives of the arguments
   expr_t composeDerivatives(expr_t darg1, expr_t darg2);
@@ -1029,6 +1017,8 @@ public:
                              map<pair<int, int>, temporary_terms_t> &temp_terms_map,
                              map<expr_t, pair<int, pair<int, int>>> &reference_count,
                              bool is_matlab) const override;
+  void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                  map<expr_t, pair<int, int>> &reference_count) const override;
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
@@ -1042,21 +1032,14 @@ public:
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
   void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                     bool lhs_rhs,
+                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                      deriv_node_temp_terms_t &tef_terms) const override;
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
-  void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const override;
   static double eval_opcode(double v1, BinaryOpcode op_code, double v2, int derivOrder) noexcept(false);
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
@@ -1152,7 +1135,7 @@ public:
 private:
   expr_t computeDerivative(int deriv_id) override;
   int cost(int cost, bool is_matlab) const override;
-  int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const override;
+  int cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const override;
   int cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const override;
   //! Returns the derivative of this node if darg1, darg2 and darg3 are the derivatives of the arguments
   expr_t composeDerivatives(expr_t darg1, expr_t darg2, expr_t darg3);
@@ -1165,6 +1148,8 @@ public:
                              map<pair<int, int>, temporary_terms_t> &temp_terms_map,
                              map<expr_t, pair<int, pair<int, int>>> &reference_count,
                              bool is_matlab) const override;
+  void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                  map<expr_t, pair<int, int>> &reference_count) const override;
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
@@ -1178,21 +1163,14 @@ public:
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
   void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                     bool lhs_rhs,
+                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                      deriv_node_temp_terms_t &tef_terms) const override;
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
-  void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const override;
   static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) noexcept(false);
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
@@ -1282,6 +1260,8 @@ public:
                              map<pair<int, int>, temporary_terms_t> &temp_terms_map,
                              map<expr_t, pair<int, pair<int, int>>> &reference_count,
                              bool is_matlab) const override;
+  void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                  map<expr_t, pair<int, int>> &reference_count) const override;
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override = 0;
   void writeJsonAST(ostream &output) const override = 0;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic = true) const override = 0;
@@ -1295,25 +1275,18 @@ public:
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic = true) const override = 0;
   void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                     bool lhs_rhs,
+                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                      deriv_node_temp_terms_t &tef_terms) const override = 0;
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override = 0;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
-  void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const override;
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
   unsigned int compileExternalFunctionArguments(ostream &CompileCode, unsigned int &instruction_number,
-                                                bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                bool lhs_rhs,
+                                                const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                                 const deriv_node_temp_terms_t &tef_terms) const;
 
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override = 0;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override = 0;
   expr_t toStatic(DataTree &static_datatree) const override = 0;
   void computeXrefs(EquationInfo &ei) const override = 0;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
@@ -1392,16 +1365,10 @@ public:
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
   void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                     bool lhs_rhs,
+                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                      deriv_node_temp_terms_t &tef_terms) const override;
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const override;
@@ -1421,18 +1388,12 @@ public:
                                  int top_level_symb_id_arg,
                                  const vector<expr_t> &arguments_arg,
                                  int inputIndex_arg);
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
   void compile(ostream &CompileCode, unsigned int &instruction_number,
-               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-               const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+               bool lhs_rhs,
+               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                const deriv_node_temp_terms_t &tef_terms) const override;
   void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                    const temporary_terms_t &temporary_terms,
@@ -1443,8 +1404,8 @@ public:
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
   void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                     bool lhs_rhs,
+                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                      deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
@@ -1467,18 +1428,12 @@ public:
                                   const vector<expr_t> &arguments_arg,
                                   int inputIndex1_arg,
                                   int inputIndex2_arg);
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
   void compile(ostream &CompileCode, unsigned int &instruction_number,
-               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-               const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+               bool lhs_rhs,
+               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                const deriv_node_temp_terms_t &tef_terms) const override;
   void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                    const temporary_terms_t &temporary_terms,
@@ -1489,8 +1444,8 @@ public:
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
   void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                     bool lhs_rhs,
+                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                      deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
@@ -1507,13 +1462,9 @@ public:
                              map<pair<int, int>, temporary_terms_t> &temp_terms_map,
                              map<expr_t, pair<int, pair<int, int>>> &reference_count,
                              bool is_matlab) const override;
+  void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                  map<expr_t, pair<int, int>> &reference_count) const override;
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   expr_t clone(DataTree &datatree) const override;
   int maxEndoLead() const override;
@@ -1551,10 +1502,9 @@ public:
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
   BinaryOpNode *normalizeEquationHelper(const set<expr_t> &contain_var, expr_t rhs) const override;
   void compile(ostream &CompileCode, unsigned int &instruction_number,
-               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-               const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+               bool lhs_rhs,
+               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                const deriv_node_temp_terms_t &tef_terms) const override;
-  void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   bool containsEndogenous() const override;
@@ -1589,13 +1539,9 @@ public:
                              map<pair<int, int>, temporary_terms_t> &temp_terms_map,
                              map<expr_t, pair<int, pair<int, int>>> &reference_count,
                              bool is_matlab) const override;
+  void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
+                                  map<expr_t, pair<int, int>> &reference_count) const override;
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
-  void computeTemporaryTerms(map<expr_t, int> &reference_count,
-                             temporary_terms_t &temporary_terms,
-                             map<expr_t, pair<int, int>> &first_occurence,
-                             int Curr_block,
-                             vector< vector<temporary_terms_t>> &v_temporary_terms,
-                             int equation) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   expr_t clone(DataTree &datatree) const override;
   int maxEndoLead() const override;
@@ -1633,10 +1579,9 @@ public:
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
   BinaryOpNode *normalizeEquationHelper(const set<expr_t> &contain_var, expr_t rhs) const override;
   void compile(ostream &CompileCode, unsigned int &instruction_number,
-               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-               const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+               bool lhs_rhs,
+               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                const deriv_node_temp_terms_t &tef_terms) const override;
-  void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   bool containsEndogenous() const override;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index f6b1b9365a2195603f8860f8460a44421036e1c9..2bd95e794082f0b86414c4db853b3a40457ac489 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -39,18 +39,6 @@ void
 ModelTree::copyHelper(const ModelTree &m)
 {
   auto f = [this](expr_t e) { return e->clone(*this); };
-  auto convert_vector_tt = [f](vector<temporary_terms_t> vtt)
-                           {
-                             vector<temporary_terms_t> vtt2;
-                             for (const auto &tt : vtt)
-                               {
-                                 temporary_terms_t tt2;
-                                 for (const auto &it : tt)
-                                   tt2.insert(f(it));
-                                 vtt2.push_back(tt2);
-                               }
-                             return vtt2;
-                           };
 
   // Equations
   for (const auto &it : m.equations)
@@ -81,16 +69,12 @@ ModelTree::copyHelper(const ModelTree &m)
                                    };
 
   // Temporary terms
-  for (const auto &it : m.temporary_terms)
-    temporary_terms.insert(f(it));
   for (const auto &it : m.temporary_terms_mlv)
     temporary_terms_mlv[f(it.first)] = f(it.second);
   for (const auto &it : m.temporary_terms_derivatives)
     temporary_terms_derivatives.push_back(convert_temporary_terms_t(it));
   for (const auto &it : m.temporary_terms_idxs)
     temporary_terms_idxs[f(it.first)] = it.second;
-  for (const auto &it : m.v_temporary_terms)
-    v_temporary_terms.push_back(convert_vector_tt(it));
   for (const auto &it : m.params_derivs_temporary_terms)
     params_derivs_temporary_terms[it.first] = convert_temporary_terms_t(it.second);
   for (const auto &it : m.params_derivs_temporary_terms_idxs)
@@ -112,6 +96,11 @@ ModelTree::copyHelper(const ModelTree &m)
         v[it2.first] = f(it2.second);
       blocks_derivatives.push_back(v);
     }
+
+  for (const auto &it : m.blocks_temporary_terms)
+    blocks_temporary_terms.push_back(convert_temporary_terms_t(it));
+  for (const auto &it : m.blocks_temporary_terms_idxs)
+    blocks_temporary_terms_idxs[f(it.first)] = it.second;
 }
 
 ModelTree::ModelTree(SymbolTable &symbol_table_arg,
@@ -136,12 +125,10 @@ ModelTree::ModelTree(const ModelTree &m) :
   equation_tags{m.equation_tags},
   computed_derivs_order{m.computed_derivs_order},
   NNZDerivatives{m.NNZDerivatives},
-  v_temporary_terms_inuse{m.v_temporary_terms_inuse},
   eq_idx_block2orig{m.eq_idx_block2orig},
   endo_idx_block2orig{m.endo_idx_block2orig},
   eq_idx_orig2block{m.eq_idx_orig2block},
   endo_idx_orig2block{m.endo_idx_orig2block},
-  map_idx{m.map_idx},
   blocks{m.blocks},
   endo2block{m.endo2block},
   eq2block{m.eq2block},
@@ -168,11 +155,8 @@ ModelTree::operator=(const ModelTree &m)
   derivatives.clear();
   params_derivatives.clear();
 
-  temporary_terms.clear();
   temporary_terms_mlv.clear();
   temporary_terms_derivatives.clear();
-  v_temporary_terms.clear();
-  v_temporary_terms_inuse = m.v_temporary_terms_inuse;
   params_derivs_temporary_terms.clear();
   params_derivs_temporary_terms_idxs.clear();
 
@@ -183,12 +167,13 @@ ModelTree::operator=(const ModelTree &m)
   endo_idx_block2orig = m.endo_idx_block2orig;
   eq_idx_orig2block = m.eq_idx_orig2block;
   endo_idx_orig2block = m.endo_idx_orig2block;
-  map_idx = m.map_idx;
   equation_type_and_normalized_equation.clear();
   blocks_derivatives.clear();
   blocks = m.blocks;
   endo2block = m.endo2block;
   eq2block = m.eq2block;
+  blocks_temporary_terms.clear();
+  blocks_temporary_terms_idxs.clear();
   is_equation_linear = m.is_equation_linear;
   endo2eq = m.endo2eq;
   cutoff = m.cutoff;
@@ -342,7 +327,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context) const
           catch (ExprNode::EvalException &e)
             {
               cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " (line " << equations_lineno[eq] << ") and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
-              d1->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModelSparse, temporary_terms, {});
+              d1->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModelSparse, {}, {});
               cerr << endl;
               exit(EXIT_FAILURE);
             }
@@ -1050,12 +1035,7 @@ ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
         else
           ++it2;
 
-  // Fill the (now obsolete) temporary_terms structure
-  temporary_terms.clear();
-  for (const auto &it : temp_terms_map)
-    temporary_terms.insert(it.second.begin(), it.second.end());
-
-  // Fill the new structure
+  // Fill the structures
   temporary_terms_derivatives.clear();
   temporary_terms_derivatives.resize(derivatives.size());
   for (int order = 0; order < static_cast<int>(derivatives.size()); order++)
@@ -1070,6 +1050,43 @@ ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
       temporary_terms_idxs[it] = idx++;
 }
 
+void
+ModelTree::computeBlockTemporaryTerms()
+{
+  int nb_blocks = blocks.size();
+  blocks_temporary_terms.resize(nb_blocks);
+
+  map<expr_t, pair<int, int>> reference_count;
+  for (int blk = 0; blk < nb_blocks; blk++)
+    {
+      for (int eq = 0; eq < blocks[blk].size; eq++)
+        {
+          if (eq < blocks[blk].getRecursiveSize() && isBlockEquationRenormalized(blk, eq))
+            getBlockEquationRenormalizedExpr(blk, eq)->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+          else
+            getBlockEquationExpr(blk, eq)->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+        }
+      for (const auto &[ignore, d] : blocks_derivatives[blk])
+        d->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+
+      additionalBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+    }
+
+  // Compute indices in the temporary terms vector
+  int idx = 0;
+  blocks_temporary_terms_idxs.clear();
+  for (int blk = 0; blk < nb_blocks; blk++)
+    for (auto tt : blocks_temporary_terms[blk])
+      blocks_temporary_terms_idxs[tt] = idx++;
+}
+
+void
+ModelTree::additionalBlockTemporaryTerms(int blk,
+                                         vector<temporary_terms_t> &blocks_temporary_terms,
+                                         map<expr_t, pair<int, int>> &reference_count) const
+{
+}
+
 void
 ModelTree::writeModelLocalVariableTemporaryTerms(temporary_terms_t &temp_term_union,
                                                  const temporary_terms_idxs_t &tt_idxs,
@@ -1297,34 +1314,28 @@ ModelTree::testNestedParenthesis(const string &str) const
 }
 
 void
-ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const
+ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, bool dynamic, bool steady_dynamic) const
 {
-  // Local var used to keep track of temp nodes already written
-  temporary_terms_t tt2;
   // To store the functions that have already been written in the form TEF* = ext_fun();
   deriv_node_temp_terms_t tef_terms;
-  for (auto it : tt)
+  for (auto [tt, idx] : temporary_terms_idxs)
     {
-      if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-        {
-          it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
-        }
+      if (dynamic_cast<AbstractExternalFunctionNode *>(tt))
+        tt->compileExternalFunctionOutput(code_file, instruction_number, false, blocks_temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 
-      FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(map_idx.find(it->idx)->second));
+      FNUMEXPR_ fnumexpr(TemporaryTerm, idx);
       fnumexpr.write(code_file, instruction_number);
-      it->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
+      tt->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
       if (dynamic)
         {
-          FSTPT_ fstpt(static_cast<int>(map_idx.find(it->idx)->second));
+          FSTPT_ fstpt(idx);
           fstpt.write(code_file, instruction_number);
         }
       else
         {
-          FSTPST_ fstpst(static_cast<int>(map_idx.find(it->idx)->second));
+          FSTPST_ fstpst(idx);
           fstpst.write(code_file, instruction_number);
         }
-      // Insert current node into tt2
-      tt2.insert(it);
     }
 }
 
@@ -1444,7 +1455,7 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type,
 }
 
 void
-ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_number, bool dynamic, bool steady_dynamic) const
 {
   for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
     {
@@ -1465,8 +1476,8 @@ ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_n
 
       if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
         {
-          lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
-          rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
+          lhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, dynamic, steady_dynamic);
+          rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, dynamic, steady_dynamic);
 
           FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
           fbinary.write(code_file, instruction_number);
@@ -1476,7 +1487,7 @@ ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_n
         }
       else // The right hand side of the equation is empty ==> residual=lhs;
         {
-          lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
+          lhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, dynamic, steady_dynamic);
           FSTPR_ fstpr(eq);
           fstpr.write(code_file, instruction_number);
         }
@@ -1876,19 +1887,19 @@ ModelTree::writeJsonModelEquations(ostream &output, bool residuals) const
         {
           output << R"({"residual": {)"
                  << R"("lhs": ")";
-          lhs->writeJsonOutput(output, temporary_terms, {});
+          lhs->writeJsonOutput(output, {}, {});
           output << R"(")";
 
           output << R"(, "rhs": ")";
-          rhs->writeJsonOutput(output, temporary_terms, {});
+          rhs->writeJsonOutput(output, {}, {});
           output << R"(")";
           try
             {
               // Test if the right hand side of the equation is empty.
-              if (rhs->eval(eval_context_t()) != 0)
+              if (rhs->eval({}) != 0)
                 {
                   output << R"(, "rhs": ")";
-                  rhs->writeJsonOutput(output, temporary_terms, {});
+                  rhs->writeJsonOutput(output, {}, {});
                   output << R"(")";
                 }
             }
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index e1700cd9137169a4a898267fd9411f1df2113e07..47a3a122578d05a70d61fb10ab6a6e10aa5d0f64 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -121,9 +121,6 @@ protected:
   then the IDs of parameters (in non-decreasing order)*/
   map<pair<int,int>, map<vector<int>, expr_t>> params_derivatives;
 
-  //! Storage for temporary terms in block/bytecode mode
-  temporary_terms_t temporary_terms;
-
   //! Used model local variables, that will be treated as temporary terms
   /*! See the comments in ModelTree::computeTemporaryTerms() */
   map<expr_t, expr_t, ExprNodeLess> temporary_terms_mlv;
@@ -135,10 +132,6 @@ protected:
   //! Stores, for each temporary term, its index in the MATLAB/Julia vector
   temporary_terms_idxs_t temporary_terms_idxs;
 
-  //! Temporary terms for block decomposed models
-  vector<vector<temporary_terms_t>> v_temporary_terms;
-  vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
-
   //! Temporary terms for parameter derivatives, under a disaggregated form
   /*! The pair of integers is to be interpreted as in param_derivatives */
   map<pair<int, int>, temporary_terms_t> params_derivs_temporary_terms;
@@ -165,8 +158,6 @@ protected:
      Set by updateReverseVariableEquationOrderings() */
   vector<int> eq_idx_orig2block, endo_idx_orig2block;
 
-  map_idx_t map_idx;
-
   //! 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;
 
@@ -202,6 +193,13 @@ protected:
      It verifies: ∀i, eq2block[endo2eq[i]] = endo2block[i] */
   vector<int> eq2block;
 
+  //! Temporary terms for block decomposed models (one block per element of the vector)
+  vector<temporary_terms_t> blocks_temporary_terms;
+
+  /* Stores, for each temporary term in block decomposed models, its index in
+     the vector of all temporary terms */
+  temporary_terms_idxs_t blocks_temporary_terms_idxs;
+
   //! the file containing the model and the derivatives code
   ofstream code_file;
 
@@ -218,13 +216,22 @@ protected:
   void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
   //! Computes temporary terms (for all equations and derivatives)
   void computeTemporaryTerms(bool is_matlab, bool no_tmp_terms);
+  //! Computes temporary terms per block
+  void computeBlockTemporaryTerms();
+  /* Add additional temporary terms for a given block. This method is called by
+     computeBlockTemporaryTerms(). It does nothing by default, but is meant to
+     be overriden by subclasses (actually by DynamicModel, who needs extra
+     temporary terms for derivatives w.r.t. exogenous and other endogenous) */
+  virtual void additionalBlockTemporaryTerms(int blk,
+                                             vector<temporary_terms_t> &blocks_temporary_terms,
+                                             map<expr_t, pair<int, int>> &reference_count) const;
   //! Computes temporary terms for the file containing parameters derivatives
   void computeParamsDerivativesTemporaryTerms();
   //! Writes temporary terms
   void writeTemporaryTerms(const temporary_terms_t &tt, temporary_terms_t &temp_term_union, const temporary_terms_idxs_t &tt_idxs, ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
   void writeJsonTemporaryTerms(const temporary_terms_t &tt, temporary_terms_t &temp_term_union, ostream &output, deriv_node_temp_terms_t &tef_terms, const string &concat) const;
   //! Compiles temporary terms
-  void compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const;
+  void compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, bool dynamic, bool steady_dynamic) const;
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File(const string &filename, int &u_count_int, bool &file_open, bool is_two_boundaries, int block_mfs) const;
   //! Fixes output when there are more than 32 nested parens, Issue #1201
@@ -245,7 +252,7 @@ protected:
   void writeJsonModelEquations(ostream &output, bool residuals) const;
   void writeJsonModelLocalVariables(ostream &output, deriv_node_temp_terms_t &tef_terms) const;
   //! Compiles model equations
-  void compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+  void compileModelEquations(ostream &code_file, unsigned int &instruction_number, bool dynamic, bool steady_dynamic) const;
 
   //! Writes LaTeX model file
   void writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const;
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 370c4bd684141eb5d18d7d4a48bf2b6159747b67..db96d48c40789087d6bafb8829433d4dbe03394b 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -29,24 +29,6 @@
 void
 StaticModel::copyHelper(const StaticModel &m)
 {
-  auto f = [this](expr_t e) { return e->clone(*this); };
-
-  auto convert_vector_tt = [f](vector<temporary_terms_t> vtt)
-                           {
-                             vector<temporary_terms_t> vtt2;
-                             for (const auto &tt : vtt)
-                               {
-                                 temporary_terms_t tt2;
-                                 for (const auto &it : tt)
-                                   tt2.insert(f(it));
-                                 vtt2.push_back(tt2);
-                               }
-                             return vtt2;
-                           };
-
-  for (const auto &it : m.v_temporary_terms_local)
-    v_temporary_terms_local.push_back(convert_vector_tt(it));
-
 }
 
 StaticModel::StaticModel(SymbolTable &symbol_table_arg,
@@ -57,9 +39,7 @@ StaticModel::StaticModel(SymbolTable &symbol_table_arg,
 }
 
 StaticModel::StaticModel(const StaticModel &m) :
-  ModelTree{m},
-  map_idx2{m.map_idx2},
-  global_temporary_terms{m.global_temporary_terms}
+  ModelTree{m}
 {
   copyHelper(m);
 }
@@ -69,10 +49,6 @@ StaticModel::operator=(const StaticModel &m)
 {
   ModelTree::operator=(m);
 
-  v_temporary_terms_local.clear();
-  map_idx2 = m.map_idx2;
-  global_temporary_terms = m.global_temporary_terms;
-
   copyHelper(m);
 
   return *this;
@@ -126,11 +102,11 @@ StaticModel::StaticModel(const DynamicModel &m) :
 }
 
 void
-StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
+StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id) const
 {
   if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) });
       it != derivatives[1].end())
-    it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
+    it->second->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
   else
     {
       FLDZ_ fldz;
@@ -139,11 +115,11 @@ StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_nu
 }
 
 void
-StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
+StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag) const
 {
   if (auto it = blocks_derivatives[blk].find({ eq, var, lag });
       it != blocks_derivatives[blk].end())
-    it->second->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
+    it->second->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
   else
     {
       FLDZ_ fldz;
@@ -151,115 +127,10 @@ StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instr
     }
 }
 
-void
-StaticModel::computeTemporaryTermsOrdered()
-{
-  map<expr_t, pair<int, int>> first_occurence;
-  map<expr_t, int> reference_count;
-  BinaryOpNode *eq_node;
-  ostringstream tmp_s;
-  map_idx.clear();
-
-  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);
-
-  v_temporary_terms_inuse = vector<temporary_terms_inuse_t>(nb_blocks);
-
-  map_idx2 = vector<map_idx_t>(nb_blocks);
-
-  temporary_terms.clear();
-
-  //local temporay terms
-  for (int block = 0; block < nb_blocks; block++)
-    {
-      map<expr_t, int> reference_count_local;
-      map<expr_t, pair<int, int>> first_occurence_local;
-      temporary_terms_t temporary_terms_l;
-
-      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++)
-        {
-          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-            getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, i);
-          else
-            {
-              eq_node = static_cast<BinaryOpNode *>(getBlockEquationExpr(block, i));
-              eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, i);
-            }
-        }
-      for (const auto &[ignore, id] : blocks_derivatives[block])
-        id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, block_size-1);
-      v_temporary_terms_inuse[block] = {};
-      computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]);
-    }
-
-  // global temporay terms
-  for (int block = 0; block < nb_blocks; block++)
-    {
-      // Compute the temporary terms reordered
-      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++)
-        {
-          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-            getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
-          else
-            {
-              eq_node = getBlockEquationExpr(block, i);
-              eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
-            }
-        }
-      for (const auto &[ignore, id] : blocks_derivatives[block])
-        id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
-    }
-
-  for (int block = 0; block < nb_blocks; block++)
-    {
-      // Collecte the temporary terms reordered
-      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++)
-        {
-          if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-            getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          else
-            {
-              eq_node = getBlockEquationExpr(block, i);
-              eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-            }
-        }
-      for (const auto &[ignore, id] : blocks_derivatives[block])
-        id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-      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;
-    }
-  computeTemporaryTermsMapping(temporary_terms, map_idx);
-}
-
-void
-StaticModel::computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, map_idx_t &map_idx)
-{
-  // Add a mapping form node ID to temporary terms order
-  int j = 0;
-  for (auto temporary_term : temporary_terms)
-    map_idx[temporary_term->idx] = j++;
-}
-
 void
 StaticModel::writeModelEquationsOrdered_M(const string &basename) const
 {
-  temporary_terms_t local_temporary_terms;
-  if (global_temporary_terms)
-    local_temporary_terms = temporary_terms;
-
+  temporary_terms_t temporary_terms; // Temp terms written so far
   constexpr ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabStaticModelSparse;
 
   //----------------------------------------------------------------------
@@ -299,40 +170,44 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
           && simulation_type != BlockSimulationType::evaluateForward)
         output << " g1 = spalloc("  << block_mfs << ", " << block_mfs << ", " << blocks_derivatives[block].size() << ");" << endl;
 
-      if (v_temporary_terms_inuse[block].size())
-        {
-          output << "  global";
-          for (int it : v_temporary_terms_inuse[block])
-            output << " T" << it;
-          output << ";" << endl;
-        }
+      // Declare global temp terms from this block and the previous ones
+      bool global_keyword_written = false;
+      for (int blk2 = 0; blk2 <= block; blk2++)
+        for (auto tt : blocks_temporary_terms[blk2])
+          {
+            if (!global_keyword_written)
+              {
+                output << "  global";
+                global_keyword_written = true;
+              }
+            output << " ";
+            tt->writeOutput(output, local_output_type, blocks_temporary_terms[blk2], {});
+          }
+      if (global_keyword_written)
+        output << ";" << endl;
 
       if (simulation_type != BlockSimulationType::evaluateBackward
           && simulation_type != BlockSimulationType::evaluateForward)
         output << "  residual=zeros(" << block_mfs << ",1);" << endl;
 
       // The equations
+
       deriv_node_temp_terms_t tef_terms;
-      for (int i = 0; i < block_size; i++)
+      for (auto it : blocks_temporary_terms[block])
         {
-          if (!global_temporary_terms)
-            local_temporary_terms = v_temporary_terms[block][i];
-          temporary_terms_t tt2;
-          if (v_temporary_terms[block].size())
-            for (auto it : v_temporary_terms[block][i])
-              {
-                if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-                  it->writeExternalFunctionOutput(output, local_output_type, tt2, {}, tef_terms);
-
-                output << "  ";
-                it->writeOutput(output, local_output_type, local_temporary_terms, {}, tef_terms);
-                output << " = ";
-                it->writeOutput(output, local_output_type, tt2, {}, tef_terms);
-                // Insert current node into tt2
-                tt2.insert(it);
-                output << ";" << endl;
-              }
+          if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+            it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, {}, tef_terms);
+
+          output << "  ";
+          it->writeOutput(output, local_output_type, blocks_temporary_terms[block], {}, tef_terms);
+          output << " = ";
+          it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms);
+          temporary_terms.insert(it);
+          output << ";" << endl;
+        }
 
+      for (int i = 0; i < block_size; i++)
+        {
           int variable_ID = getBlockVariableID(block, i);
           int equation_ID = getBlockEquationID(block, i);
           EquationType equ_type = getBlockEquationType(block, i);
@@ -340,7 +215,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
           BinaryOpNode *eq_node = getBlockEquationExpr(block, i);
           expr_t lhs = eq_node->arg1, rhs = eq_node->arg2;
           ostringstream tmp_output;
-          lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms, {});
+          lhs->writeOutput(tmp_output, local_output_type, temporary_terms, {});
           switch (simulation_type)
             {
             case BlockSimulationType::evaluateBackward:
@@ -350,16 +225,16 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
               if (equ_type == EquationType::evaluate)
                 {
                   output << tmp_output.str() << " = ";
-                  rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+                  rhs->writeOutput(output, local_output_type, temporary_terms, {});
                 }
               else if (equ_type == EquationType::evaluate_s)
                 {
                   eq_node = getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  lhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+                  lhs->writeOutput(output, local_output_type, temporary_terms, {});
                   output << " = ";
-                  rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+                  rhs->writeOutput(output, local_output_type, temporary_terms, {});
                 }
               else
                 {
@@ -376,7 +251,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
                 goto evaluation;
               output << "  " << "residual(" << i+1-block_recursive << ") = ("
                      << tmp_output.str() << ") - (";
-              rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
+              rhs->writeOutput(output, local_output_type, temporary_terms, {});
               output << ");" << endl;
               break;
             default:
@@ -395,7 +270,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
             {
               auto [eq, var, ignore] = indices;
               output << "    g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = ";
-              id->writeOutput(output, local_output_type, local_temporary_terms, {});
+              id->writeOutput(output, local_output_type, temporary_terms, {});
               output << ";" << endl;
             }
           break;
@@ -408,9 +283,8 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
 }
 
 void
-StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx) const
+StaticModel::writeModelEquationsCode(const string &basename) const
 {
-
   ostringstream tmp_output;
   ofstream code_file;
   unsigned int instruction_number = 0;
@@ -431,6 +305,11 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
   Write_Inf_To_Bin_File(basename + "/model/bytecode/static.bin", u_count_int, file_open, false, symbol_table.endo_nbr());
   file_open = true;
 
+  // Compute the union of temporary terms from residuals and 1st derivatives
+  temporary_terms_t temporary_terms = temporary_terms_derivatives[0];
+  copy(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end(),
+       inserter(temporary_terms, temporary_terms.end()));
+
   //Temporary variables declaration
   FDIMST_ fdimst(temporary_terms.size());
   fdimst.write(code_file, instruction_number);
@@ -448,13 +327,9 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
                            symbol_table.endo_nbr());
   fbeginblock.write(code_file, instruction_number);
 
-  // Add a mapping form node ID to temporary terms order
-  int j = 0;
-  for (auto temporary_term : temporary_terms)
-    map_idx[temporary_term->idx] = j++;
-  compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, false, false);
+  compileTemporaryTerms(code_file, instruction_number, false, false);
 
-  compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, false, false);
+  compileModelEquations(code_file, instruction_number, false, false);
 
   FENDEQU_ fendequ;
   fendequ.write(code_file, instruction_number);
@@ -481,7 +356,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
             my_derivatives[eq].clear();
           my_derivatives[eq].emplace_back(var, count_u);
 
-          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
+          d1->compile(code_file, instruction_number, false, temporary_terms_idxs, false, false);
 
           FSTPSU_ fstpsu(count_u);
           fstpsu.write(code_file, instruction_number);
@@ -543,7 +418,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
             my_derivatives[eq].clear();
           my_derivatives[eq].emplace_back(var, count_u);
 
-          d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
+          d1->compile(code_file, instruction_number, false, temporary_terms_idxs, false, false);
           FSTPG2_ fstpg2(eq, var);
           fstpg2.write(code_file, instruction_number);
         }
@@ -564,7 +439,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
 }
 
 void
-StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map_idx, vector<map_idx_t> map_idx2) const
+StaticModel::writeModelEquationsCode_Block(const string &basename) const
 {
   struct Uff_l
   {
@@ -587,7 +462,6 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
   Uff Uf[symbol_table.endo_nbr()];
   map<expr_t, int> reference_count;
   vector<int> feedback_variables;
-  deriv_node_temp_terms_t tef_terms;
   bool file_open = false;
 
   filesystem::create_directories(basename + "/model/bytecode");
@@ -601,7 +475,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
     }
   //Temporary variables declaration
 
-  FDIMST_ fdimst(temporary_terms.size());
+  FDIMST_ fdimst(blocks_temporary_terms_idxs.size());
   fdimst.write(code_file, instruction_number);
 
   for (int block = 0; block < static_cast<int>(blocks.size()); block++)
@@ -649,27 +523,22 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
       fjmp_if_eval.write(code_file, instruction_number);
       int prev_instruction_number = instruction_number;
 
-      for (i = 0; i < block_size; i++)
+      //The Temporary terms
+      deriv_node_temp_terms_t tef_terms;
+      for (auto it : blocks_temporary_terms[block])
         {
-          //The Temporary terms
-          temporary_terms_t tt2;
-          if (v_temporary_terms[block].size())
-            {
-              for (auto it : v_temporary_terms[block][i])
-                {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-                    it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
-
-                  FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(map_idx.find(it->idx)->second));
-                  fnumexpr.write(code_file, instruction_number);
-                  it->compile(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
-                  FSTPST_ fstpst(static_cast<int>(map_idx.find(it->idx)->second));
-                  fstpst.write(code_file, instruction_number);
-                  // Insert current node into tt2
-                  tt2.insert(it);
-                }
-            }
+          if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+            it->compileExternalFunctionOutput(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false, tef_terms);
 
+          FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+          fnumexpr.write(code_file, instruction_number);
+          it->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false, tef_terms);
+          FSTPST_ fstpst(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+          fstpst.write(code_file, instruction_number);
+        }
+
+      for (i = 0; i < block_size; i++)
+        {
           // The equations
           int variable_ID, equation_ID;
           EquationType equ_type;
@@ -688,16 +557,16 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
                   eq_node = getBlockEquationExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
+                  rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
+                  lhs->compile(code_file, instruction_number, true, blocks_temporary_terms_idxs, false, false);
                 }
               else if (equ_type == EquationType::evaluate_s)
                 {
                   eq_node = getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
+                  rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
+                  lhs->compile(code_file, instruction_number, true, blocks_temporary_terms_idxs, false, false);
                 }
               break;
             case BlockSimulationType::solveBackwardComplete:
@@ -716,8 +585,8 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
               eq_node = getBlockEquationExpr(block, i);
               lhs = eq_node->arg1;
               rhs = eq_node->arg2;
-              lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
-              rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
+              lhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
+              rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
 
               FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
               fbinary.write(code_file, instruction_number);
@@ -741,7 +610,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
                 FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
                 fnumexpr.write(code_file, instruction_number);
               }
-              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx, temporary_terms);
+              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0));
               {
                 FSTPG_ fstpg(0);
                 fstpg.write(code_file, instruction_number);
@@ -773,7 +642,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
                       Uf[eqr].Ufl->var = varr;
                       FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
                       fnumexpr.write(code_file, instruction_number);
-                      compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, map_idx, temporary_terms);
+                      compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0);
                       FSTPSU_ fstpsu(count_u);
                       fstpsu.write(code_file, instruction_number);
                       count_u++;
@@ -836,31 +705,21 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
       code_file.seekp(pos3);
       prev_instruction_number = instruction_number;
 
-      temporary_terms_t tt2, tt3;
-      deriv_node_temp_terms_t tef_terms2;
-
-      for (i = 0; i < block_size; i++)
+      tef_terms.clear();
+      for (auto it : blocks_temporary_terms[block])
         {
-          if (v_temporary_terms_local[block].size())
-            {
-              for (const auto &it : v_temporary_terms_local[block][i])
-                {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-                    it->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
+          if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+            it->compileExternalFunctionOutput(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false, tef_terms);
 
-                  FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(map_idx2[block].find(it->idx)->second));
-                  fnumexpr.write(code_file, instruction_number);
-
-                  it->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);
-
-                  FSTPST_ fstpst(static_cast<int>(map_idx2[block].find(it->idx)->second));
-                  fstpst.write(code_file, instruction_number);
-                  // Insert current node into tt2
-                  tt3.insert(it);
-                  tt2.insert(it);
-                }
-            }
+          FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+          fnumexpr.write(code_file, instruction_number);
+          it->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false, tef_terms);
+          FSTPST_ fstpst(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+          fstpst.write(code_file, instruction_number);
+        }
 
+      for (i = 0; i < block_size; i++)
+        {
           // The equations
           int variable_ID, equation_ID;
           EquationType equ_type;
@@ -879,16 +738,16 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
                   eq_node = getBlockEquationExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
-                  lhs->compile(code_file, instruction_number, true, tt2, map_idx2[block], false, false);
+                  rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
+                  lhs->compile(code_file, instruction_number, true, blocks_temporary_terms_idxs, false, false);
                 }
               else if (equ_type == EquationType::evaluate_s)
                 {
                   eq_node = getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
-                  lhs->compile(code_file, instruction_number, true, tt2, map_idx2[block], false, false);
+                  rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
+                  lhs->compile(code_file, instruction_number, true, blocks_temporary_terms_idxs, false, false);
                 }
               break;
             case BlockSimulationType::solveBackwardComplete:
@@ -907,8 +766,8 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
               eq_node = getBlockEquationExpr(block, i);
               lhs = eq_node->arg1;
               rhs = eq_node->arg2;
-              lhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
-              rhs->compile(code_file, instruction_number, false, tt2, map_idx2[block], false, false);
+              lhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
+              rhs->compile(code_file, instruction_number, false, blocks_temporary_terms_idxs, false, false);
 
               FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
               fbinary.write(code_file, instruction_number);
@@ -929,7 +788,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
             FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
             fnumexpr.write(code_file, instruction_number);
           }
-          compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx2[block], tt2 /*temporary_terms*/);
+          compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0));
           {
             FSTPG2_ fstpg2(0, 0);
             fstpg2.write(code_file, instruction_number);
@@ -948,7 +807,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
               FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0);
               fnumexpr.write(code_file, instruction_number);
 
-              compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, map_idx2[block], tt2 /*temporary_terms*/);
+              compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0);
 
               FSTPG2_ fstpg2(eq, var);
               fstpg2.write(code_file, instruction_number);
@@ -1086,15 +945,12 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
 
       determineLinearBlocks();
 
-      global_temporary_terms = true;
       if (!no_tmp_terms)
-        computeTemporaryTermsOrdered();
+        computeBlockTemporaryTerms();
     }
   else
     {
       computeTemporaryTerms(true, no_tmp_terms);
-      if (bytecode && !no_tmp_terms)
-        computeTemporaryTermsMapping(temporary_terms, map_idx);
 
       /* Must be called after computeTemporaryTerms(), because it depends on
          temporary_terms_mlv to be filled */
@@ -1854,9 +1710,9 @@ void
 StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const
 {
   if (block && bytecode)
-    writeModelEquationsCode_Block(basename, map_idx, map_idx2);
+    writeModelEquationsCode_Block(basename);
   else if (!block && bytecode)
-    writeModelEquationsCode(basename, map_idx);
+    writeModelEquationsCode(basename);
   else if (block && !bytecode)
     {
       writeModelEquationsOrdered_M(basename);
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index a4683c75395fd1585a2a0f96c3520bce8bc85066..1c191feb4b42da5308e808daf90be43b347eae71 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -33,9 +33,6 @@ class DynamicModel;
 class StaticModel : public ModelTree
 {
 private:
-  //! local temporary terms for block decomposed models
-  vector<vector<temporary_terms_t>> v_temporary_terms_local;
-
   //! Writes static model file (standard Matlab version)
   void writeStaticMFile(const string &basename) const;
 
@@ -55,10 +52,10 @@ private:
   void writeModelEquationsOrdered_M(const string &basename) const;
 
   //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
-  void writeModelEquationsCode_Block(const string &basename, map_idx_t map_idx, vector<map_idx_t> map_idx2) const;
+  void writeModelEquationsCode_Block(const string &basename) const;
 
   //! Writes the code of the model in virtual machine bytecode
-  void writeModelEquationsCode(const string &basename, map_idx_t map_idx) const;
+  void writeModelEquationsCode(const string &basename) const;
 
   //! Computes jacobian and prepares for equation normalization
   /*! Using values from initval/endval blocks and parameter initializations:
@@ -67,17 +64,10 @@ private:
   */
   void evaluateJacobian(const eval_context_t &eval_context, jacob_map_t *j_m, bool dynamic);
 
-  vector<map_idx_t> map_idx2;
-
-  //! sorts the temporary terms in the blocks order
-  void computeTemporaryTermsOrdered();
-  //! creates a mapping from the index of temporary terms to a natural index
-  void computeTemporaryTermsMapping(temporary_terms_t &temporary_terms, map_idx_t &map_idx);
-
   //! Write derivative code of an equation w.r. to a variable
-  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
+  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id) const;
   //! Write chain rule derivative code of an equation w.r. to a variable
-  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const;
+  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag) const;
 
   //! Get the type corresponding to a derivation ID
   SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
@@ -90,9 +80,6 @@ private:
   //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
   void computeChainRuleJacobian();
 
-  //! Indicate if the temporary terms are computed for the overall model (true) or not (false). Default value true
-  bool global_temporary_terms{true};
-
   //! Helper functions for writeStaticModel
   void writeStaticModelHelper(const string &basename,
                               const string &name, const string &retvalname,