diff --git a/src/Bytecode.hh b/src/Bytecode.hh
index c7bef351aad70cf35f758c801a11395ec7caa34b..9d867e55c95065cebfad731e13a3757bb06e0143 100644
--- a/src/Bytecode.hh
+++ b/src/Bytecode.hh
@@ -928,7 +928,7 @@ public:
   }
   FBEGINBLOCK_(int size_arg, BlockSimulationType type_arg, int first_element, int block_size,
                const vector<int> &variable_arg, const vector<int> &equation_arg,
-               bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg, int nb_col_jacob_arg,
+               bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int u_count_int_arg, int nb_col_jacob_arg,
                int det_exo_size_arg, int nb_col_det_exo_jacob_arg, int exo_size_arg, int nb_col_exo_jacob_arg, int other_endo_size_arg, int nb_col_other_endo_jacob_arg,
                vector<int> det_exogenous_arg, vector<int> exogenous_arg, vector<int> other_endogenous_arg) :
     BytecodeInstruction{Tags::FBEGINBLOCK},
@@ -955,7 +955,7 @@ public:
   }
   FBEGINBLOCK_(int size_arg, BlockSimulationType type_arg, int first_element, int block_size,
                const vector<int> &variable_arg, const vector<int> &equation_arg,
-               bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg, int nb_col_jacob_arg) :
+               bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int u_count_int_arg, int nb_col_jacob_arg) :
     BytecodeInstruction{Tags::FBEGINBLOCK},
     size{size_arg},
     type{type_arg},
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 8174bc13d3425c4348c8483f7b717bd8e156cd96..a7b140dcbf1c5f894a0c0f5781786eecb0daba9b 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -165,26 +165,6 @@ DynamicModel::operator=(const DynamicModel &m)
   return *this;
 }
 
-void
-DynamicModel::writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
-{
-  if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) });
-      it != derivatives[1].end())
-    it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms, temporary_terms_idxs, tef_terms);
-  else
-    code_file << FLDZ_{};
-}
-
-void
-DynamicModel::writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
-{
-  if (auto it = blocks_derivatives[blk].find({ eq, var, lag });
-      it != blocks_derivatives[blk].end())
-    it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms, temporary_terms_idxs, tef_terms);
-  else
-    code_file << FLDZ_{};
-}
-
 void
 DynamicModel::additionalBlockTemporaryTerms(int blk,
                                             vector<vector<temporary_terms_t>> &blocks_temporary_terms,
@@ -332,6 +312,39 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
     }
 }
 
+void
+DynamicModel::writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block,
+                                                      const temporary_terms_t &temporary_terms_union,
+                                                      const deriv_node_temp_terms_t &tef_terms) const
+{
+  constexpr ExprNodeBytecodeOutputType output_type {ExprNodeBytecodeOutputType::dynamicModel};
+
+  /* FIXME: there is an inconsistency between endos and the following 3 other
+     variable types. For the latter, the index of equation within the block is
+     taken from FNUMEXPR, while it is taken from FSTPG3 for the former. */
+  for (const auto &[indices, d] : blocks_derivatives_exo[block])
+    {
+      const auto &[eq, var, lag] {indices};
+      code_file << FNUMEXPR_{ExpressionType::FirstExoDerivative, eq, 0, lag};
+      d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+      code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo[block].at({ var, lag })};
+    }
+  for (const auto &[indices, d] : blocks_derivatives_exo_det[block])
+    {
+      const auto &[eq, var, lag] {indices};
+      code_file << FNUMEXPR_{ExpressionType::FirstExodetDerivative, eq, 0, lag};
+      d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+      code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo_det[block].at({ var, lag })};
+    }
+  for (const auto &[indices, d] : blocks_derivatives_other_endo[block])
+    {
+      const auto &[eq, var, lag] {indices};
+      code_file << FNUMEXPR_{ExpressionType::FirstOtherEndoDerivative, eq, 0, lag};
+      d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+      code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_other_endo[block].at({ var, lag })};
+    }
+}
+
 void
 DynamicModel::writeDynamicPerBlockCFiles(const string &basename) const
 {
@@ -591,301 +604,56 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
 void
 DynamicModel::writeDynamicBlockBytecode(const string &basename) const
 {
-  struct Uff_l
-  {
-    int u, var, lag;
-    Uff_l *pNext;
-  };
-
-  struct Uff
-  {
-    Uff_l *Ufl, *Ufl_First;
-  };
-
-  int i, v;
-  string tmp_s;
-  ostringstream tmp_output;
-  expr_t lhs = nullptr, rhs = nullptr;
-  BinaryOpNode *eq_node;
-  Uff Uf[symbol_table.endo_nbr()];
-  map<expr_t, int> reference_count;
-  vector<int> feedback_variables;
-  bool file_open = false;
-  BytecodeWriter code_file{basename + "/model/bytecode/dynamic.cod"};
+  BytecodeWriter code_file {basename + "/model/bytecode/dynamic.cod"};
 
-  //Temporary variables declaration
+  const string bin_filename {basename + "/model/bytecode/dynamic.bin"};
+  ofstream bin_file {bin_filename, ios::out | ios::binary};
+  if (!bin_file.is_open())
+    {
+      cerr << R"(Error : Can't open file ")" << bin_filename << R"(" for writing)" << endl;
+      exit(EXIT_FAILURE);
+    }
 
+  // Temporary variables declaration
   code_file << FDIMT_{static_cast<int>(blocks_temporary_terms_idxs.size())};
 
-  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
+  for (int block {0}; block < static_cast<int>(blocks.size()); block++)
     {
-      feedback_variables.clear();
-      if (block > 0)
-        code_file << FENDBLOCK_{};
-      int count_u;
-      int u_count_int = 0;
-      BlockSimulationType simulation_type = blocks[block].simulation_type;
-      int block_size = blocks[block].size;
-      int block_mfs = blocks[block].mfs_size;
-      int block_recursive = blocks[block].getRecursiveSize();
-      int block_max_lag = blocks[block].max_lag;
-      int block_max_lead = blocks[block].max_lead;
-
-      if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple
-          || simulation_type == BlockSimulationType::solveTwoBoundariesComplete
-          || simulation_type == BlockSimulationType::solveBackwardComplete
-          || simulation_type == BlockSimulationType::solveForwardComplete)
-        {
-          writeBlockBytecodeBinFile(basename, block, u_count_int, file_open,
-                                    simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple);
-          file_open = true;
-        }
-
-      code_file << FBEGINBLOCK_{block_mfs,
-          simulation_type,
-          blocks[block].first_equation,
-          block_size,
-          endo_idx_block2orig,
-          eq_idx_block2orig,
-          blocks[block].linear,
-          symbol_table.endo_nbr(),
-          block_max_lag,
-          block_max_lead,
-          u_count_int,
-          static_cast<int>(blocks_jacob_cols_endo[block].size()),
-          static_cast<int>(blocks_exo_det[block].size()),
-          static_cast<int>(blocks_jacob_cols_exo_det[block].size()),
-          static_cast<int>(blocks_exo[block].size()),
-          static_cast<int>(blocks_jacob_cols_exo[block].size()),
-          static_cast<int>(blocks_other_endo[block].size()),
-          static_cast<int>(blocks_jacob_cols_other_endo[block].size()),
-          vector<int>(blocks_exo_det[block].begin(), blocks_exo_det[block].end()),
-          vector<int>(blocks_exo[block].begin(), blocks_exo[block].end()),
-          vector<int>(blocks_other_endo[block].begin(), blocks_other_endo[block].end())};
-
-      temporary_terms_t temporary_terms_union;
-
-      //The Temporary terms
-      deriv_node_temp_terms_t tef_terms;
-
-      auto write_eq_tt = [&](int eq)
-                         {
-                           for (auto it : blocks_temporary_terms[block][eq])
-                             {
-                               if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-                                 it->writeBytecodeExternalFunctionOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-
-                               code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)};
-                               it->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                               code_file << FSTPT_{blocks_temporary_terms_idxs.at(it)};
-                               temporary_terms_union.insert(it);
-                             }
-                         };
-
-      // The equations
-      for (i = 0; i < block_size; i++)
-        {
-          write_eq_tt(i);
-
-          int variable_ID, equation_ID;
-          EquationType equ_type;
-
-          switch (simulation_type)
-            {
-            evaluation:
-            case BlockSimulationType::evaluateBackward:
-            case BlockSimulationType::evaluateForward:
-              equ_type = getBlockEquationType(block, i);
-              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
-              if (equ_type == EquationType::evaluate)
-                {
-                  eq_node = getBlockEquationExpr(block, i);
-                  lhs = eq_node->arg1;
-                  rhs = eq_node->arg2;
-                  rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                  lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                }
-              else if (equ_type == EquationType::evaluateRenormalized)
-                {
-                  eq_node = getBlockEquationRenormalizedExpr(block, i);
-                  lhs = eq_node->arg1;
-                  rhs = eq_node->arg2;
-                  rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                  lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                }
-              break;
-            case BlockSimulationType::solveBackwardComplete:
-            case BlockSimulationType::solveForwardComplete:
-            case BlockSimulationType::solveTwoBoundariesComplete:
-            case BlockSimulationType::solveTwoBoundariesSimple:
-              if (i < block_recursive)
-                goto evaluation;
-              variable_ID = getBlockVariableID(block, i);
-              equation_ID = getBlockEquationID(block, i);
-              feedback_variables.push_back(variable_ID);
-              Uf[equation_ID].Ufl = nullptr;
-              goto end;
-            default:
-            end:
-              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
-              eq_node = getBlockEquationExpr(block, i);
-              lhs = eq_node->arg1;
-              rhs = eq_node->arg2;
-              lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-              rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-
-              code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
-            }
-        }
-      code_file << FENDEQU_{};
-
-      /* If the block is not of type “evaluate backward/forward”, then we write
-         the temporary terms for derivatives at this point, i.e. before the
-         JMPIFEVAL, because they will be needed in both “simulate” and
-         “evaluate” modes. */
-      if (simulation_type != BlockSimulationType::evaluateBackward
-          && simulation_type != BlockSimulationType::evaluateForward)
-        write_eq_tt(blocks[block].size);
-
-      // Get the current code_file position and jump if evaluating
-      int pos_jmpifeval = code_file.getInstructionCounter();
-      code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being
-
-      /* Write the derivatives for the “simulate” mode (not needed if the block
-         is of type “evaluate backward/forward”) */
-      if (simulation_type != BlockSimulationType::evaluateBackward
-          && simulation_type != BlockSimulationType::evaluateForward)
-        {
-          switch (simulation_type)
-            {
-            case BlockSimulationType::solveBackwardSimple:
-            case BlockSimulationType::solveForwardSimple:
-              code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0};
-              writeBytecodeDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-              code_file << FSTPG_{0};
-              break;
-
-            case BlockSimulationType::solveBackwardComplete:
-            case BlockSimulationType::solveForwardComplete:
-            case BlockSimulationType::solveTwoBoundariesComplete:
-            case BlockSimulationType::solveTwoBoundariesSimple:
-              count_u = feedback_variables.size();
-              for (const auto &[indices, ignore] : blocks_derivatives[block])
-                {
-                  auto [eq, var, lag] = indices;
-                  int eqr = getBlockEquationID(block, eq);
-                  int varr = getBlockVariableID(block, var);
-                  if (eq >= block_recursive and var >= block_recursive)
-                    {
-                      if (lag != 0
-                          && (simulation_type == BlockSimulationType::solveForwardComplete
-                              || simulation_type == BlockSimulationType::solveBackwardComplete))
-                        continue;
-                      if (!Uf[eqr].Ufl)
-                        {
-                          Uf[eqr].Ufl = static_cast<Uff_l *>(malloc(sizeof(Uff_l)));
-                          Uf[eqr].Ufl_First = Uf[eqr].Ufl;
-                        }
-                      else
-                        {
-                          Uf[eqr].Ufl->pNext = static_cast<Uff_l *>(malloc(sizeof(Uff_l)));
-                          Uf[eqr].Ufl = Uf[eqr].Ufl->pNext;
-                        }
-                      Uf[eqr].Ufl->pNext = nullptr;
-                      Uf[eqr].Ufl->u = count_u;
-                      Uf[eqr].Ufl->var = varr;
-                      Uf[eqr].Ufl->lag = lag;
-                      code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag};
-                      writeBytecodeChainRuleDerivative(code_file, block, eq, var, lag, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                      code_file << FSTPU_{count_u};
-                      count_u++;
-                    }
-                }
-              for (i = 0; i < block_size; i++)
-                {
-                  if (i >= block_recursive)
-                    {
-                      code_file << FLDR_{i-block_recursive} << FLDZ_{};
-
-                      v = getBlockEquationID(block, i);
-                      for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
-                        code_file << FLDU_{Uf[v].Ufl->u}
-                          << FLDV_{SymbolType::endogenous, Uf[v].Ufl->var, Uf[v].Ufl->lag}
-                          << FBINARY_{BinaryOpcode::times}
-                          << FCUML_{};
-                      Uf[v].Ufl = Uf[v].Ufl_First;
-                      while (Uf[v].Ufl)
-                        {
-                          Uf[v].Ufl_First = Uf[v].Ufl->pNext;
-                          free(Uf[v].Ufl);
-                          Uf[v].Ufl = Uf[v].Ufl_First;
-                        }
-                      code_file << FBINARY_{BinaryOpcode::minus}
-                        << FSTPU_{i - block_recursive};
-                    }
-                }
-              break;
-            default:
-              break;
-            }
-        }
-
-      // Jump unconditionally after the block
-      int pos_jmp = code_file.getInstructionCounter();
-      code_file << FJMP_{0}; // Use 0 as jump offset for the time being
-      // Update jump offset for previous JMPIFEVAL
-      code_file.overwriteInstruction(pos_jmpifeval, FJMPIFEVAL_{pos_jmp-pos_jmpifeval});
-
-      /* If the block is of type “evaluate backward/forward”, then write the
-         temporary terms for derivatives at this point, because they have not
-         been written before the JMPIFEVAL. */
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        write_eq_tt(blocks[block].size);
-
-      // Write the derivatives for the “evaluate” mode
-      for (const auto &[indices, d] : blocks_derivatives[block])
-        {
-          auto [eq, var, lag] = indices;
-          int eqr = getBlockEquationID(block, eq);
-          int varr = getBlockVariableID(block, var);
-          code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag};
-          writeBytecodeDerivative(code_file, eqr, varr, lag, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-          code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_endo[block].at({ var, lag })};
-        }
-      for (const auto &[indices, d] : blocks_derivatives_exo[block])
-        {
-          auto [eqr, var, lag] = indices;
-          int eq = getBlockEquationID(block, eqr);
-          int varr = 0; // Dummy value, actually unused by the bytecode MEX
-          code_file << FNUMEXPR_{ExpressionType::FirstExoDerivative, eqr, varr, lag};
-          d->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-          code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo[block].at({ var, lag })};
-        }
-      for (const auto &[indices, d] : blocks_derivatives_exo_det[block])
-        {
-          auto [eqr, var, lag] = indices;
-          int eq = getBlockEquationID(block, eqr);
-          int varr = 0; // Dummy value, actually unused by the bytecode MEX
-          code_file << FNUMEXPR_{ExpressionType::FirstExodetDerivative, eqr, varr, lag};
-          d->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-          code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo_det[block].at({ var, lag })};
-        }
-      for (const auto &[indices, d] : blocks_derivatives_other_endo[block])
-        {
-          auto [eqr, var, lag] = indices;
-          int eq = getBlockEquationID(block, eqr);
-          int varr = 0; // Dummy value, actually unused by the bytecode MEX
-          code_file << FNUMEXPR_{ExpressionType::FirstOtherEndoDerivative, eqr, varr, lag};
-          d->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-          code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_other_endo[block].at({ var, lag })};
-        }
-
-      // Update jump offset for previous JMP
-      int pos_end_block = code_file.getInstructionCounter();
-      code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1});
+      const BlockSimulationType simulation_type {blocks[block].simulation_type};
+
+      // Write section of .bin file except for evaluate blocks and solve simple blocks
+      const int u_count {simulation_type == BlockSimulationType::solveTwoBoundariesSimple
+                         || simulation_type == BlockSimulationType::solveTwoBoundariesComplete
+                         || simulation_type == BlockSimulationType::solveBackwardComplete
+                         || simulation_type == BlockSimulationType::solveForwardComplete
+                         ? writeBlockBytecodeBinFile(bin_file, block)
+                         : 0};
+
+      code_file << FBEGINBLOCK_{blocks[block].mfs_size,
+                                simulation_type,
+                                blocks[block].first_equation,
+                                blocks[block].size,
+                                endo_idx_block2orig,
+                                eq_idx_block2orig,
+                                blocks[block].linear,
+                                symbol_table.endo_nbr(),
+                                blocks[block].max_lag,
+                                blocks[block].max_lead,
+                                u_count,
+                                static_cast<int>(blocks_jacob_cols_endo[block].size()),
+                                static_cast<int>(blocks_exo_det[block].size()),
+                                static_cast<int>(blocks_jacob_cols_exo_det[block].size()),
+                                static_cast<int>(blocks_exo[block].size()),
+                                static_cast<int>(blocks_jacob_cols_exo[block].size()),
+                                static_cast<int>(blocks_other_endo[block].size()),
+                                static_cast<int>(blocks_jacob_cols_other_endo[block].size()),
+                                { blocks_exo_det[block].begin(), blocks_exo_det[block].end() },
+                                { blocks_exo[block].begin(), blocks_exo[block].end() },
+                                { blocks_other_endo[block].begin(), blocks_other_endo[block].end() }};
+
+      writeBlockBytecodeHelper<true>(code_file, block);
     }
-  code_file << FENDBLOCK_{} << FEND_{};
+  code_file << FEND_{};
 }
 
 void
@@ -1333,60 +1101,6 @@ DynamicModel::printNonZeroHessianEquations(ostream &output) const
     output << "]";
 }
 
-void
-DynamicModel::writeBlockBytecodeBinFile(const string &basename, int num, int &u_count_int,
-                                        bool &file_open, bool is_two_boundaries) const
-{
-  int j;
-  std::ofstream SaveCode;
-  string filename = basename + "/model/bytecode/dynamic.bin";
-
-  if (file_open)
-    SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate);
-  else
-    SaveCode.open(filename, ios::out | ios::binary);
-  if (!SaveCode.is_open())
-    {
-      cerr << R"(Error : Can't open file ")" << filename << R"(" for writing)" << endl;
-      exit(EXIT_FAILURE);
-    }
-  u_count_int = 0;
-  int block_size = blocks[num].size;
-  int block_mfs = blocks[num].mfs_size;
-  int block_recursive = blocks[num].getRecursiveSize();
-  for (const auto &[indices, ignore] : blocks_derivatives[num])
-    {
-      auto [eq, var, lag] = indices;
-      if (lag != 0 && !is_two_boundaries)
-        continue;
-      if (eq >= block_recursive && var >= block_recursive)
-        {
-          int v = eq - block_recursive;
-          SaveCode.write(reinterpret_cast<char *>(&v), sizeof(v));
-          int varr = var - block_recursive + lag * block_mfs;
-          SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
-          SaveCode.write(reinterpret_cast<const char *>(&lag), sizeof(lag));
-          int u = u_count_int + block_mfs;
-          SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
-          u_count_int++;
-        }
-    }
-
-  if (is_two_boundaries)
-    u_count_int += block_mfs;
-  for (j = block_recursive; j < block_size; j++)
-    {
-      int varr = getBlockVariableID(num, j);
-      SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
-    }
-  for (j = block_recursive; j < block_size; j++)
-    {
-      int eqr = getBlockEquationID(num, j);
-      SaveCode.write(reinterpret_cast<char *>(&eqr), sizeof(eqr));
-    }
-  SaveCode.close();
-}
-
 void
 DynamicModel::writeDynamicBlockMFile(const string &basename) const
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 4d44f8836002a2e2292974a9916ea625d4fc1b74..491cbeef8011d62db3a36ace7d786fbb720d53ab 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -140,11 +140,12 @@ private:
   void writeDynamicPerBlockCFiles(const string &basename) const;
   //! Writes the code of the block-decomposed model in virtual machine bytecode
   void writeDynamicBlockBytecode(const string &basename) const;
+  // Writes derivatives w.r.t. exo, exo det and other endogenous
+  void writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block,
+                                               const temporary_terms_t &temporary_terms_union,
+                                               const deriv_node_temp_terms_t &tef_terms) const override;
   //! Writes the code of the model in virtual machine bytecode
   void writeDynamicBytecode(const string &basename) const;
-  //! Adds per-block information for bytecode simulation in a separate .bin file
-  void writeBlockBytecodeBinFile(const string &basename, int num, int &u_count_int, bool &file_open,
-                                 bool is_two_boundaries) const;
 
   void writeSetAuxiliaryVariables(const string &basename, bool julia) const;
   void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const;
@@ -175,11 +176,6 @@ private:
                                      vector<vector<temporary_terms_t>> &blocks_temporary_terms,
                                      map<expr_t, tuple<int, int, int>> &reference_count) const override;
 
-  //! Write derivative bytecode of an equation w.r. to a variable
-  void writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
-  //! Write chain rule derivative bytecode of an equation w.r. to a variable
-  void writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
-
   SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
   int getLagByDerivID(int deriv_id) const noexcept(false) override;
   int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;
@@ -295,6 +291,12 @@ private:
      (with a lag), and also as a contemporaneous diff lag auxvar). */
   vector<int> getVARDerivIDs(int lhs_symb_id, int lead_lag) const;
 
+  int
+  getBlockJacobianEndoCol(int blk, int var, int lag) const override
+  {
+    return blocks_jacob_cols_endo[blk].at({ var, lag });
+  }
+
 public:
   DynamicModel(SymbolTable &symbol_table_arg,
                NumericalConstants &num_constants_arg,
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 0e70d8dae79d671dcd40df3043dd87513dbd206c..b03581cae97366298dd502a9ad398faa6872ea52 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1262,6 +1262,49 @@ ModelTree::writeBytecodeBinFile(const string &filename, bool is_two_boundaries)
   return u_count;
 }
 
+int
+ModelTree::writeBlockBytecodeBinFile(ofstream &bin_file, int block) const
+{
+  int u_count {0};
+  int block_size {blocks[block].size};
+  int block_mfs {blocks[block].mfs_size};
+  int block_recursive {blocks[block].getRecursiveSize()};
+  BlockSimulationType simulation_type {blocks[block].simulation_type};
+  bool is_two_boundaries {simulation_type == BlockSimulationType::solveTwoBoundariesComplete
+                          || simulation_type == BlockSimulationType::solveTwoBoundariesSimple};
+  for (const auto &[indices, ignore] : blocks_derivatives[block])
+    {
+      const auto &[eq, var, lag] {indices};
+      if (lag != 0 && !is_two_boundaries)
+        continue;
+      if (eq >= block_recursive && var >= block_recursive)
+        {
+          int v {eq - block_recursive};
+          bin_file.write(reinterpret_cast<char *>(&v), sizeof v);
+          int varr {var - block_recursive + lag * block_mfs};
+          bin_file.write(reinterpret_cast<char *>(&varr), sizeof varr);
+          bin_file.write(reinterpret_cast<const char *>(&lag), sizeof lag);
+          int u {u_count + block_mfs};
+          bin_file.write(reinterpret_cast<char *>(&u), sizeof u);
+          u_count++;
+        }
+    }
+
+  if (is_two_boundaries)
+    u_count += block_mfs;
+  for (int j {block_recursive}; j < block_size; j++)
+    {
+      int varr {getBlockVariableID(block, j)};
+      bin_file.write(reinterpret_cast<char *>(&varr), sizeof varr);
+    }
+  for (int j {block_recursive}; j < block_size; j++)
+    {
+      int eqr {getBlockEquationID(block, j)};
+      bin_file.write(reinterpret_cast<char *>(&eqr), sizeof eqr);
+    }
+  return u_count;
+}
+
 void
 ModelTree::writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const
 {
@@ -1823,3 +1866,11 @@ ModelTree::getRHSFromLHS(expr_t lhs) const
       return eq->arg2;
   throw ExprNode::MatchFailureException{"Cannot find an equation with the requested LHS"};
 }
+
+void
+ModelTree::writeBlockBytecodeAdditionalDerivatives([[maybe_unused]] BytecodeWriter &code_file,
+                                                   [[maybe_unused]] int block,
+                                                   [[maybe_unused]] const temporary_terms_t &temporary_terms_union,
+                                                   [[maybe_unused]] const deriv_node_temp_terms_t &tef_terms) const
+{
+}
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index a195cf79093cfd790322aae67ebf3cd706d588c6..f11ac7c9ab8b655a9085928f6dc79cf51350cdbc 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -251,6 +251,9 @@ protected:
      file.
      Returns the number of first derivatives w.r.t. endogenous variables */
   int writeBytecodeBinFile(const string &filename, bool is_two_boundaries) const;
+  //! Adds per-block information for bytecode simulation in a separate .bin file
+  int writeBlockBytecodeBinFile(ofstream &bin_file, int block) const;
+
   //! Fixes output when there are more than 32 nested parens, Issue #1201
   void fixNestedParenthesis(ostringstream &output, map<string, string> &tmp_paren_vars, bool &message_printed) const;
   //! Tests if string contains more than 32 nested parens, Issue #1201
@@ -279,6 +282,17 @@ protected:
   template<bool dynamic>
   void writeBytecodeHelper(BytecodeWriter &code_file) const;
 
+  // Helper for writing blocks in bytecode
+  template<bool dynamic>
+  void writeBlockBytecodeHelper(BytecodeWriter &code_file, int block) const;
+
+  /* Write additional derivatives w.r.t. to exogenous, exogenous det and other endo
+     in block+bytecode mode. Does nothing by default, but overriden by
+     DynamicModel which needs those. */
+  virtual void writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block,
+                                                       const temporary_terms_t &temporary_terms_union,
+                                                       const deriv_node_temp_terms_t &tef_terms) const;
+
   /* Helper for writing JSON output for residuals and derivatives.
      Returns mlv and derivatives output at each derivation order. */
   template<bool dynamic>
@@ -430,6 +444,10 @@ protected:
       Assumes that derivatives have already been computed. */
   map<tuple<int, int, int>, expr_t> collectFirstOrderDerivativesEndogenous();
 
+  /* Get column number within Jacobian of a given block.
+     “var” is the block-specific endogenous variable index. */
+  virtual int getBlockJacobianEndoCol(int blk, int var, int lag) const = 0;
+
 private:
   //! Internal helper for the copy constructor and assignment operator
   /*! Copies all the structures that contain ExprNode*, by the converting the
@@ -1205,6 +1223,226 @@ ModelTree::writeBytecodeHelper(BytecodeWriter &code_file) const
   code_file << FENDBLOCK_{} << FEND_{};
 }
 
+template<bool dynamic>
+void
+ModelTree::writeBlockBytecodeHelper(BytecodeWriter &code_file, int block) const
+{
+  constexpr ExprNodeBytecodeOutputType output_type
+    { dynamic ? ExprNodeBytecodeOutputType::dynamicModel : ExprNodeBytecodeOutputType::staticModel };
+  constexpr ExprNodeBytecodeOutputType assignment_lhs_output_type
+    { dynamic ? ExprNodeBytecodeOutputType::dynamicAssignmentLHS : ExprNodeBytecodeOutputType::staticAssignmentLHS };
+
+  const BlockSimulationType simulation_type {blocks[block].simulation_type};
+  const int block_size {blocks[block].size};
+  const int block_mfs {blocks[block].mfs_size};
+  const int block_recursive {blocks[block].getRecursiveSize()};
+
+  temporary_terms_t temporary_terms_union;
+  deriv_node_temp_terms_t tef_terms;
+
+  auto write_eq_tt = [&](int eq)
+  {
+    for (auto it : blocks_temporary_terms[block][eq])
+      {
+        if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+          it->writeBytecodeExternalFunctionOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+
+        code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)};
+        it->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+        if constexpr(dynamic)
+          code_file << FSTPT_{blocks_temporary_terms_idxs.at(it)};
+        else
+          code_file << FSTPST_{blocks_temporary_terms_idxs.at(it)};
+        temporary_terms_union.insert(it);
+      }
+  };
+
+  // The equations
+  for (int i {0}; i < block_size; i++)
+    {
+      write_eq_tt(i);
+
+      switch (simulation_type)
+        {
+        evaluation:
+        case BlockSimulationType::evaluateBackward:
+        case BlockSimulationType::evaluateForward:
+          code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
+          if (EquationType equ_type {getBlockEquationType(block, i)};
+              equ_type == EquationType::evaluate)
+            {
+              BinaryOpNode *eq_node {getBlockEquationExpr(block, i)};
+              expr_t lhs {eq_node->arg1};
+              expr_t rhs {eq_node->arg2};
+              rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+              lhs->writeBytecodeOutput(code_file, assignment_lhs_output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+            }
+          else if (equ_type == EquationType::evaluateRenormalized)
+            {
+              BinaryOpNode *eq_node {getBlockEquationRenormalizedExpr(block, i)};
+              expr_t lhs {eq_node->arg1};
+              expr_t rhs {eq_node->arg2};
+              rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+              lhs->writeBytecodeOutput(code_file, assignment_lhs_output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+            }
+          break;
+        case BlockSimulationType::solveBackwardComplete:
+        case BlockSimulationType::solveForwardComplete:
+        case BlockSimulationType::solveTwoBoundariesComplete:
+        case BlockSimulationType::solveTwoBoundariesSimple:
+          if (i < block_recursive)
+            goto evaluation;
+          [[fallthrough]];
+        default:
+          code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
+          BinaryOpNode *eq_node {getBlockEquationExpr(block, i)};
+          expr_t lhs {eq_node->arg1};
+          expr_t rhs {eq_node->arg2};
+          lhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+          rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+          code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
+        }
+    }
+  code_file << FENDEQU_{};
+
+  /* If the block is not of type “evaluate backward/forward”, then we write
+     the temporary terms for derivatives at this point, i.e. before the
+     JMPIFEVAL, because they will be needed in both “simulate” and
+     “evaluate” modes. */
+  if (simulation_type != BlockSimulationType::evaluateBackward
+      && simulation_type != BlockSimulationType::evaluateForward)
+    write_eq_tt(blocks[block].size);
+
+  // Get the current code_file position and jump if evaluating
+  int pos_jmpifeval {code_file.getInstructionCounter()};
+  code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being
+
+  /* Write the derivatives for the “simulate” mode (not needed if the block
+     is of type “evaluate backward/forward”) */
+  if (simulation_type != BlockSimulationType::evaluateBackward
+      && simulation_type != BlockSimulationType::evaluateForward)
+    {
+      switch (simulation_type)
+        {
+        case BlockSimulationType::solveBackwardSimple:
+        case BlockSimulationType::solveForwardSimple:
+          {
+            int eqr {getBlockEquationID(block, 0)};
+            int varr {getBlockVariableID(block, 0)};
+            code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, 0};
+            // Get contemporaneous derivative of the single variable in the block
+            if (auto it { blocks_derivatives[block].find({ 0, 0, 0 }) };
+                it != blocks_derivatives[block].end())
+              it->second->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+            else
+              code_file << FLDZ_{};
+            code_file << FSTPG_{0};
+          }
+          break;
+
+        case BlockSimulationType::solveBackwardComplete:
+        case BlockSimulationType::solveForwardComplete:
+        case BlockSimulationType::solveTwoBoundariesComplete:
+        case BlockSimulationType::solveTwoBoundariesSimple:
+          {
+            // For each equation, stores a list of tuples (index_u, var, lag)
+            vector<vector<tuple<int, int, int>>> Uf(symbol_table.endo_nbr());
+
+            for (int count_u {block_mfs};
+                 const auto &[indices, ignore] : blocks_derivatives[block])
+              {
+                const auto &[eq, var, lag] {indices};
+                int eqr {getBlockEquationID(block, eq)};
+                int varr {getBlockVariableID(block, var)};
+                if (eq >= block_recursive && var >= block_recursive)
+                  {
+                    if constexpr(dynamic)
+                      if (lag != 0
+                          && (simulation_type == BlockSimulationType::solveForwardComplete
+                              || simulation_type == BlockSimulationType::solveBackwardComplete))
+                        continue;
+                    code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag};
+                    if (auto it { blocks_derivatives[block].find({ eq, var, lag }) };
+                        it != blocks_derivatives[block].end())
+                      it->second->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+                    else
+                      code_file << FLDZ_{};
+                    if constexpr(dynamic)
+                      code_file << FSTPU_{count_u};
+                    else
+                      code_file << FSTPSU_{count_u};
+                    Uf[eqr].emplace_back(count_u, varr, lag);
+                    count_u++;
+                  }
+              }
+            for (int i {0}; i < block_size; i++)
+              if (i >= block_recursive)
+                {
+                  code_file << FLDR_{i-block_recursive} << FLDZ_{};
+
+                  int eqr {getBlockEquationID(block, i)};
+                  for (const auto &[index_u, var, lag] : Uf[eqr])
+                    {
+                      if constexpr(dynamic)
+                        code_file << FLDU_{index_u}
+                                  << FLDV_{SymbolType::endogenous, var, lag};
+                      else
+                        code_file << FLDSU_{index_u}
+                                  << FLDSV_{SymbolType::endogenous, var};
+                      code_file << FBINARY_{BinaryOpcode::times}
+                                << FCUML_{};
+                    }
+                  code_file << FBINARY_{BinaryOpcode::minus};
+                  if constexpr(dynamic)
+                    code_file << FSTPU_{i - block_recursive};
+                  else
+                    code_file << FSTPSU_{i - block_recursive};
+                }
+          }
+          break;
+        default:
+          break;
+        }
+    }
+
+  // Jump unconditionally after the block
+  int pos_jmp {code_file.getInstructionCounter()};
+  code_file << FJMP_{0}; // Use 0 as jump offset for the time being
+  // Update jump offset for previous JMPIFEVAL
+  code_file.overwriteInstruction(pos_jmpifeval, FJMPIFEVAL_{pos_jmp-pos_jmpifeval});
+
+  /* If the block is of type “evaluate backward/forward”, then write the
+     temporary terms for derivatives at this point, because they have not
+     been written before the JMPIFEVAL. */
+  if (simulation_type == BlockSimulationType::evaluateBackward
+      || simulation_type == BlockSimulationType::evaluateForward)
+    write_eq_tt(blocks[block].size);
+
+  // Write the derivatives for the “evaluate” mode
+  for (const auto &[indices, d] : blocks_derivatives[block])
+    {
+      const auto &[eq, var, lag] {indices};
+      int eqr {getBlockEquationID(block, eq)};
+      int varr {getBlockVariableID(block, var)};
+      code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag};
+      d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
+      if constexpr(dynamic)
+        code_file << FSTPG3_{eq, var, lag, getBlockJacobianEndoCol(block, var, lag)};
+      else
+        code_file << FSTPG2_{eq, getBlockJacobianEndoCol(block, var, lag)};
+    }
+
+  /* Write derivatives w.r.t. exo, exo det and other endogenous, but only in
+     dynamic mode */
+  writeBlockBytecodeAdditionalDerivatives(code_file, block, temporary_terms_union, tef_terms);
+
+  // Update jump offset for previous JMP
+  int pos_end_block {code_file.getInstructionCounter()};
+  code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1});
+
+  code_file << FENDBLOCK_{};
+}
+
 template<bool dynamic>
 pair<ostringstream, vector<ostringstream>>
 ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index ad32e00720783e0e4724c22dd14e78a257d33074..45d8723cb0eaf73a5ecbd2f6c59b3c78443a16a6 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -94,26 +94,6 @@ StaticModel::StaticModel(const DynamicModel &m) :
   user_set_compiler = m.user_set_compiler;
 }
 
-void
-StaticModel::writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
-{
-  if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) });
-      it != derivatives[1].end())
-    it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms, temporary_terms_idxs, tef_terms);
-  else
-    code_file << FLDZ_{};
-}
-
-void
-StaticModel::writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
-{
-  if (auto it = blocks_derivatives[blk].find({ eq, var, lag });
-      it != blocks_derivatives[blk].end())
-    it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms, temporary_terms_idxs, tef_terms);
-  else
-    code_file << FLDZ_{};
-}
-
 void
 StaticModel::writeStaticPerBlockMFiles(const string &basename) const
 {
@@ -288,368 +268,45 @@ StaticModel::writeStaticBytecode(const string &basename) const
 void
 StaticModel::writeStaticBlockBytecode(const string &basename) const
 {
-  struct Uff_l
-  {
-    int u, var, lag;
-    Uff_l *pNext;
-  };
-
-  struct Uff
-  {
-    Uff_l *Ufl, *Ufl_First;
-  };
-
-  int i, v;
-  string tmp_s;
-  ostringstream tmp_output;
-  expr_t lhs = nullptr, rhs = nullptr;
-  BinaryOpNode *eq_node;
-  Uff Uf[symbol_table.endo_nbr()];
-  map<expr_t, int> reference_count;
-  vector<int> feedback_variables;
-  bool file_open = false;
-
-  BytecodeWriter code_file{basename + "/model/bytecode/static.cod"};
-
-  //Temporary variables declaration
-  code_file << FDIMST_{static_cast<int>(blocks_temporary_terms_idxs.size())};
-
-  temporary_terms_t temporary_terms_union;
-
-  for (int block = 0; block < static_cast<int>(blocks.size()); block++)
-    {
-      feedback_variables.clear();
-      if (block > 0)
-        code_file << FENDBLOCK_{};
-      int count_u;
-      int u_count_int = 0;
-      BlockSimulationType simulation_type = blocks[block].simulation_type;
-      int block_size = blocks[block].size;
-      int block_mfs = blocks[block].mfs_size;
-      int block_recursive = blocks[block].getRecursiveSize();
-
-      if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple
-          || simulation_type == BlockSimulationType::solveTwoBoundariesComplete
-          || simulation_type == BlockSimulationType::solveBackwardComplete
-          || simulation_type == BlockSimulationType::solveForwardComplete)
-        {
-          writeBlockBytecodeBinFile(basename, block, u_count_int, file_open);
-          file_open = true;
-        }
-
-      code_file << FBEGINBLOCK_{block_mfs,
-          simulation_type,
-          blocks[block].first_equation,
-          block_size,
-          endo_idx_block2orig,
-          eq_idx_block2orig,
-          blocks[block].linear,
-          symbol_table.endo_nbr(),
-          0,
-          0,
-          u_count_int,
-          block_size};
-
-      // Get the current code_file position and jump if evaluating
-      int pos_jmpifeval = code_file.getInstructionCounter();
-      code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being
-
-      //The Temporary terms
-      deriv_node_temp_terms_t tef_terms;
-      /* Keep a backup of temporary_terms_union here, since temp. terms are
-         written a second time below. This is probably unwanted… */
-      temporary_terms_t ttu_old = temporary_terms_union;
-
-      auto write_eq_tt = [&](int eq)
-                         {
-                           for (auto it : blocks_temporary_terms[block][eq])
-                             {
-                               if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-                                 it->writeBytecodeExternalFunctionOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-
-                               code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)};
-                               it->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                               code_file << FSTPST_{blocks_temporary_terms_idxs.at(it)};
-                               temporary_terms_union.insert(it);
-                             }
-                         };
-
-      for (i = 0; i < block_size; i++)
-        {
-          write_eq_tt(i);
-
-          // The equations
-          int variable_ID, equation_ID;
-          EquationType equ_type;
-          switch (simulation_type)
-            {
-            evaluation:
-            case BlockSimulationType::evaluateBackward:
-            case BlockSimulationType::evaluateForward:
-              equ_type = getBlockEquationType(block, i);
-              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
-              if (equ_type == EquationType::evaluate)
-                {
-                  eq_node = getBlockEquationExpr(block, i);
-                  lhs = eq_node->arg1;
-                  rhs = eq_node->arg2;
-                  rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                  lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                }
-              else if (equ_type == EquationType::evaluateRenormalized)
-                {
-                  eq_node = getBlockEquationRenormalizedExpr(block, i);
-                  lhs = eq_node->arg1;
-                  rhs = eq_node->arg2;
-                  rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                  lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                }
-              break;
-            case BlockSimulationType::solveBackwardComplete:
-            case BlockSimulationType::solveForwardComplete:
-              if (i < block_recursive)
-                goto evaluation;
-              variable_ID = getBlockVariableID(block, i);
-              equation_ID = getBlockEquationID(block, i);
-              feedback_variables.push_back(variable_ID);
-              Uf[equation_ID].Ufl = nullptr;
-              goto end;
-            default:
-            end:
-              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
-              eq_node = getBlockEquationExpr(block, i);
-              lhs = eq_node->arg1;
-              rhs = eq_node->arg2;
-              lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-              rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-
-              code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
-            }
-        }
-      code_file << FENDEQU_{};
-
-      // The Jacobian if we have to solve the block
-      if (simulation_type != BlockSimulationType::evaluateBackward
-          && simulation_type != BlockSimulationType::evaluateForward)
-        {
-          // Write temporary terms for derivatives
-          write_eq_tt(blocks[block].size);
-
-          switch (simulation_type)
-            {
-            case BlockSimulationType::solveBackwardSimple:
-            case BlockSimulationType::solveForwardSimple:
-              code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, 0, 0};
-              writeBytecodeDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-              code_file << FSTPG_{0};
-              break;
-
-            case BlockSimulationType::solveBackwardComplete:
-            case BlockSimulationType::solveForwardComplete:
-              count_u = feedback_variables.size();
-              for (const auto &[indices, ignore2] : blocks_derivatives[block])
-                {
-                  auto [eq, var, ignore] = indices;
-                  int eqr = getBlockEquationID(block, eq);
-                  int varr = getBlockVariableID(block, var);
-                  if (eq >= block_recursive && var >= block_recursive)
-                    {
-                      if (!Uf[eqr].Ufl)
-                        {
-                          Uf[eqr].Ufl = static_cast<Uff_l *>(malloc(sizeof(Uff_l)));
-                          Uf[eqr].Ufl_First = Uf[eqr].Ufl;
-                        }
-                      else
-                        {
-                          Uf[eqr].Ufl->pNext = static_cast<Uff_l *>(malloc(sizeof(Uff_l)));
-                          Uf[eqr].Ufl = Uf[eqr].Ufl->pNext;
-                        }
-                      Uf[eqr].Ufl->pNext = nullptr;
-                      Uf[eqr].Ufl->u = count_u;
-                      Uf[eqr].Ufl->var = varr;
-                      code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr};
-                      writeBytecodeChainRuleDerivative(code_file, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                      code_file << FSTPSU_{count_u};
-                      count_u++;
-                    }
-                }
-              for (i = 0; i < block_size; i++)
-                if (i >= block_recursive)
-                  {
-                    code_file << FLDR_{i-block_recursive} << FLDZ_{};
-
-                    v = getBlockEquationID(block, i);
-                    for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
-                      code_file << FLDSU_{Uf[v].Ufl->u}
-                        << FLDSV_{SymbolType::endogenous, Uf[v].Ufl->var}
-                        << FBINARY_{BinaryOpcode::times}
-                        << FCUML_{};
-                    Uf[v].Ufl = Uf[v].Ufl_First;
-                    while (Uf[v].Ufl)
-                      {
-                        Uf[v].Ufl_First = Uf[v].Ufl->pNext;
-                        free(Uf[v].Ufl);
-                        Uf[v].Ufl = Uf[v].Ufl_First;
-                      }
-                    code_file << FBINARY_{BinaryOpcode::minus}
-                      << FSTPSU_{i - block_recursive};
-                  }
-              break;
-            default:
-              break;
-            }
-        }
-
-      // Jump unconditionally after the block
-      int pos_jmp = code_file.getInstructionCounter();
-      code_file << FJMP_{0}; // Use 0 as jump offset for the time being
-      // Update jump offset for previous JMPIFEVAL
-      code_file.overwriteInstruction(pos_jmpifeval, FJMPIFEVAL_{pos_jmp-pos_jmpifeval});
-
-      tef_terms.clear();
-      temporary_terms_union = ttu_old;
-
-      for (i = 0; i < block_size; i++)
-        {
-          write_eq_tt(i);
-
-          // The equations
-          int variable_ID, equation_ID;
-          EquationType equ_type;
-          switch (simulation_type)
-            {
-            evaluation_l:
-            case BlockSimulationType::evaluateBackward:
-            case BlockSimulationType::evaluateForward:
-              equ_type = getBlockEquationType(block, i);
-              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
-              if (equ_type == EquationType::evaluate)
-                {
-                  eq_node = getBlockEquationExpr(block, i);
-                  lhs = eq_node->arg1;
-                  rhs = eq_node->arg2;
-                  rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                  lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                }
-              else if (equ_type == EquationType::evaluateRenormalized)
-                {
-                  eq_node = getBlockEquationRenormalizedExpr(block, i);
-                  lhs = eq_node->arg1;
-                  rhs = eq_node->arg2;
-                  rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                  lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                }
-              break;
-            case BlockSimulationType::solveBackwardComplete:
-            case BlockSimulationType::solveForwardComplete:
-              if (i < block_recursive)
-                goto evaluation_l;
-              variable_ID = getBlockVariableID(block, i);
-              equation_ID = getBlockEquationID(block, i);
-              feedback_variables.push_back(variable_ID);
-              Uf[equation_ID].Ufl = nullptr;
-              goto end_l;
-            default:
-            end_l:
-              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
-              eq_node = getBlockEquationExpr(block, i);
-              lhs = eq_node->arg1;
-              rhs = eq_node->arg2;
-              lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-              rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-
-              code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
-            }
-        }
-      code_file << FENDEQU_{};
-
-      // The Jacobian if we have to solve the block determinsitic bloc
-
-      // Write temporary terms for derivatives
-      write_eq_tt(blocks[block].size);
-
-      switch (simulation_type)
-        {
-        case BlockSimulationType::solveBackwardSimple:
-        case BlockSimulationType::solveForwardSimple:
-          code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, 0, 0};
-          writeBytecodeDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-          code_file << FSTPG2_{0, 0};
-          break;
-        case BlockSimulationType::evaluateBackward:
-        case BlockSimulationType::evaluateForward:
-        case BlockSimulationType::solveBackwardComplete:
-        case BlockSimulationType::solveForwardComplete:
-          count_u = feedback_variables.size();
-          for (const auto &[indices, ignore2] : blocks_derivatives[block])
-            {
-              auto &[eq, var, ignore] = indices;
-              int eqr = getBlockEquationID(block, eq);
-              int varr = getBlockVariableID(block, var);
-              code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, 0};
-              writeBytecodeChainRuleDerivative(code_file, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-              code_file << FSTPG2_{eq, var};
-            }
-          break;
-        default:
-          break;
-        }
-      // Set codefile position to previous JMP_ and set the number of instructions to jump
-      // Update jump offset for previous JMP
-      int pos_end_block = code_file.getInstructionCounter();
-      code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1});
-    }
-  code_file << FENDBLOCK_{} << FEND_{};
-}
+  BytecodeWriter code_file {basename + "/model/bytecode/static.cod"};
 
-void
-StaticModel::writeBlockBytecodeBinFile(const string &basename, int num,
-                                       int &u_count_int, bool &file_open) const
-{
-  int j;
-  std::ofstream SaveCode;
-  string filename = basename + "/model/bytecode/static.bin";
-  if (file_open)
-    SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate);
-  else
-    SaveCode.open(filename, ios::out | ios::binary);
-  if (!SaveCode.is_open())
+  const string bin_filename {basename + "/model/bytecode/static.bin"};
+  ofstream bin_file {bin_filename, ios::out | ios::binary};
+  if (!bin_file.is_open())
     {
-      cerr << "Error : Can't open file " << filename << " for writing" << endl;
+      cerr << R"(Error : Can't open file ")" << bin_filename << R"(" for writing)" << endl;
       exit(EXIT_FAILURE);
     }
-  u_count_int = 0;
-  int block_size = blocks[num].size;
-  int block_mfs = blocks[num].mfs_size;
-  int block_recursive = blocks[num].getRecursiveSize();
-  for (const auto &[indices, ignore2] : blocks_derivatives[num])
-    {
-      auto [eq, var, ignore] = indices;
-      int lag = 0;
-      if (eq >= block_recursive && var >= block_recursive)
-        {
-          int v = eq - block_recursive;
-          SaveCode.write(reinterpret_cast<char *>(&v), sizeof(v));
-          int varr = var - block_recursive;
-          SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
-          SaveCode.write(reinterpret_cast<char *>(&lag), sizeof(lag));
-          int u = u_count_int + block_mfs;
-          SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
-          u_count_int++;
-        }
-    }
 
-  for (j = block_recursive; j < block_size; j++)
-    {
-      int varr = getBlockVariableID(num, j);
-      SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
-    }
-  for (j = block_recursive; j < block_size; j++)
+  // Temporary variables declaration
+  code_file << FDIMST_{static_cast<int>(blocks_temporary_terms_idxs.size())};
+
+  for (int block {0}; block < static_cast<int>(blocks.size()); block++)
     {
-      int eqr = getBlockEquationID(num, j);
-      SaveCode.write(reinterpret_cast<char *>(&eqr), sizeof(eqr));
+      const BlockSimulationType simulation_type {blocks[block].simulation_type};
+      const int block_size {blocks[block].size};
+
+      const int u_count {simulation_type == BlockSimulationType::solveBackwardComplete
+                         || simulation_type == BlockSimulationType::solveForwardComplete
+                         ? writeBlockBytecodeBinFile(bin_file, block)
+                         : 0};
+
+      code_file << FBEGINBLOCK_{blocks[block].mfs_size,
+                                simulation_type,
+                                blocks[block].first_equation,
+                                block_size,
+                                endo_idx_block2orig,
+                                eq_idx_block2orig,
+                                blocks[block].linear,
+                                symbol_table.endo_nbr(),
+                                0,
+                                0,
+                                u_count,
+                                block_size};
+
+      writeBlockBytecodeHelper<false>(code_file, block);
     }
-  SaveCode.close();
+  code_file << FEND_{};
 }
 
 void
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 8948d846d0d2756528901479a7bd410ff8900d24..facb4e09d2d53e7be5de34ee48ae53fa3308bca6 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -65,10 +65,6 @@ private:
   //! Writes the code of the model in virtual machine bytecode
   void writeStaticBytecode(const string &basename) const;
 
-  //! Adds per-block information for bytecode simulation in a separate .bin file
-  void writeBlockBytecodeBinFile(const string &basename, int num,
-                                 int &u_count_int, bool &file_open) const;
-
   //! Computes jacobian and prepares for equation normalization
   /*! Using values from initval/endval blocks and parameter initializations:
     - computes the jacobian for the model w.r. to contemporaneous variables
@@ -76,11 +72,6 @@ private:
   */
   void evaluateJacobian(const eval_context_t &eval_context, jacob_map_t *j_m, bool dynamic);
 
-  //! Write derivative bytecode of an equation w.r. to a variable
-  void writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
-  //! Write chain rule derivative bytecode of an equation w.r. to a variable
-  void writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
-
   SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
   int getLagByDerivID(int deriv_id) const noexcept(false) override;
   int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;
@@ -113,6 +104,12 @@ private:
   //! Create a legacy *_static.m file for MATLAB/Octave not yet using the temporary terms array interface
   void writeStaticMCompatFile(const string &name) const;
 
+  int
+  getBlockJacobianEndoCol([[maybe_unused]] int blk, int var, [[maybe_unused]] int lag) const override
+  {
+    return var;
+  }
+
 public:
   StaticModel(SymbolTable &symbol_table_arg,
               NumericalConstants &num_constants,