diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 61295f9065d9fc37f3992d033482d18b4a05d203..5a0d951f258eca83ea46da2aafeacc40698d46bf 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -531,202 +531,61 @@ DynamicModel::writeDynamicPerBlockCFiles(const string &basename) const
 void
 DynamicModel::writeDynamicBytecode(const string &basename) const
 {
-  ostringstream tmp_output;
-  BytecodeWriter code_file{basename + "/model/bytecode/dynamic.cod"};
-  bool file_open = false;
-
-  int count_u;
-  int u_count_int = 0;
+  // Determine the type of model (used for typing the single block)
   BlockSimulationType simulation_type;
-  if ((max_endo_lag > 0) && (max_endo_lead > 0))
+  if (max_endo_lag > 0 && max_endo_lead > 0)
     simulation_type = BlockSimulationType::solveTwoBoundariesComplete;
-  else if ((max_endo_lag >= 0) && (max_endo_lead == 0))
+  else if (max_endo_lag >= 0 && max_endo_lead == 0)
     simulation_type = BlockSimulationType::solveForwardComplete;
   else
     simulation_type = BlockSimulationType::solveBackwardComplete;
 
-  writeBytecodeBinFile(basename + "/model/bytecode/dynamic.bin", u_count_int, file_open, simulation_type == BlockSimulationType::solveTwoBoundariesComplete);
-  file_open = true;
-
-  //Temporary variables declaration
-  code_file << FDIMT_{static_cast<int>(temporary_terms_idxs.size())};
+  // First write the .bin file
+  int u_count_int { writeBytecodeBinFile(basename + "/model/bytecode/dynamic.bin",
+                                         simulation_type == BlockSimulationType::solveTwoBoundariesComplete) };
 
-  vector<int> exo, exo_det, other_endo;
+  BytecodeWriter code_file {basename + "/model/bytecode/dynamic.cod"};
 
-  for (int i = 0; i < symbol_table.exo_det_nbr(); i++)
-    exo_det.push_back(i);
-  for (int i = 0; i < symbol_table.exo_nbr(); i++)
-    exo.push_back(i);
+  // Declare temporary terms
+  code_file << FDIMT_{static_cast<int>(temporary_terms_derivatives[0].size()
+                                       + temporary_terms_derivatives[1].size())};
 
-  map<tuple<int, int, int>, expr_t> first_derivatives_reordered_endo;
-  map<tuple<int, SymbolType, int, int>, expr_t> first_derivatives_reordered_exo;
-  for (const auto & [indices, d1] : derivatives[1])
-    {
-      int deriv_id = indices[1];
-      int eq = indices[0];
-      int var { getTypeSpecificIDByDerivID(deriv_id) };
-      int lag = getLagByDerivID(deriv_id);
-      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
-        first_derivatives_reordered_endo[{ lag, var, eq }] = d1;
-      else if (getTypeByDerivID(deriv_id) == SymbolType::exogenous || getTypeByDerivID(deriv_id) == SymbolType::exogenousDet)
-        first_derivatives_reordered_exo[{ lag, getTypeByDerivID(deriv_id), var, eq }] = d1;
-    }
-  int prev_var = -1;
-  int prev_lag = -999999999;
-  int count_col_endo = 0;
-  for (const auto &it : first_derivatives_reordered_endo)
-    {
-      int var, lag;
-      tie(lag, var, ignore) = it.first;
-      if (prev_var != var || prev_lag != lag)
-        {
-          prev_var = var;
-          prev_lag = lag;
-          count_col_endo++;
-        }
-    }
-  prev_var = -1;
-  prev_lag = -999999999;
-  SymbolType prev_type{SymbolType::unusedEndogenous}; // Any non-exogenous type would do here
-  int count_col_exo = 0;
-  int count_col_det_exo = 0;
+  // Declare the (single) block
+  vector<int> exo(symbol_table.exo_nbr()), exo_det(symbol_table.exo_det_nbr());
+  iota(exo.begin(), exo.end(), 0);
+  iota(exo_det.begin(), exo_det.end(), 0);
 
-  for (const auto &it : first_derivatives_reordered_exo)
-    {
-      int var, lag;
-      SymbolType type;
-      tie(lag, type, var, ignore) = it.first;
-      if (prev_var != var || prev_lag != lag || prev_type != type)
-        {
-          prev_var = var;
-          prev_lag = lag;
-          prev_type = type;
-          if (type == SymbolType::exogenous)
-            count_col_exo++;
-          else if (type == SymbolType::exogenousDet)
-            count_col_det_exo++;
-        }
-    }
+  int jacobian_ncols_endo
+    { static_cast<int>(count_if(dyn_jacobian_cols_table.begin(), dyn_jacobian_cols_table.end(),
+                                [this](const auto &v)
+                                { return getTypeByDerivID(v.first) == SymbolType::endogenous; }))
+    };
+  int jacobian_ncols_exo {symbol_table.exo_nbr()};
+  int jacobian_ncols_exo_det {symbol_table.exo_det_nbr()};
 
   code_file << FBEGINBLOCK_{symbol_table.endo_nbr(),
-      simulation_type,
-      0,
-      symbol_table.endo_nbr(),
-      endo_idx_block2orig,
-      eq_idx_block2orig,
-      false,
-      symbol_table.endo_nbr(),
-      max_endo_lag,
-      max_endo_lead,
-      u_count_int,
-      count_col_endo,
-      symbol_table.exo_det_nbr(),
-      count_col_det_exo,
-      symbol_table.exo_nbr(),
-      count_col_exo,
-      0,
-      0,
-      exo_det,
-      exo,
-      other_endo};
-
-  temporary_terms_t temporary_terms_union;
-  deriv_node_temp_terms_t tef_terms;
-
-  writeBytecodeTemporaryTerms(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-
-  writeBytecodeModelEquations(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-
-  code_file << FENDEQU_{};
-
-  // 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
-
-  vector<vector<tuple<int, int, int>>> my_derivatives(symbol_table.endo_nbr());;
-  count_u = symbol_table.endo_nbr();
-  for (const auto & [indices, d1] : derivatives[1])
-    {
-      int deriv_id = indices[1];
-      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
-        {
-          int eq = indices[0];
-          int var { getTypeSpecificIDByDerivID(deriv_id) };
-          int lag = getLagByDerivID(deriv_id);
-          code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, var, lag};
-          if (!my_derivatives[eq].size())
-            my_derivatives[eq].clear();
-          my_derivatives[eq].emplace_back(var, lag, count_u);
-          d1->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-
-          code_file << FSTPU_{count_u};
-          count_u++;
-        }
-    }
-  for (int i = 0; i < symbol_table.endo_nbr(); i++)
-    {
-      code_file << FLDR_{i};
-      if (my_derivatives[i].size())
-        {
-          for (auto it = my_derivatives[i].begin(); it != my_derivatives[i].end(); ++it)
-            {
-              code_file << FLDU_{get<2>(*it)}
-                << FLDV_{SymbolType::endogenous, get<0>(*it), get<1>(*it)}
-                << FBINARY_{BinaryOpcode::times};
-              if (it != my_derivatives[i].begin())
-                code_file << FBINARY_{BinaryOpcode::plus};
-            }
-          code_file << FBINARY_{BinaryOpcode::minus};
-        }
-      code_file << FSTPU_{i};
-    }
-
-  // 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});
-
-  // The Jacobian
-  prev_var = -1;
-  prev_lag = -999999999;
-  count_col_endo = 0;
-  for (const auto &it : first_derivatives_reordered_endo)
-    {
-      auto [lag, var, eq] = it.first;
-      expr_t d1 = it.second;
-      code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, var, lag};
-      if (prev_var != var || prev_lag != lag)
-        {
-          prev_var = var;
-          prev_lag = lag;
-          count_col_endo++;
-        }
-      d1->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-      code_file << FSTPG3_{eq, var, lag, count_col_endo-1};
-    }
-  prev_var = -1;
-  prev_lag = -999999999;
-  count_col_exo = 0;
-  for (const auto &it : first_derivatives_reordered_exo)
-    {
-      auto [lag, ignore, var, eq] = it.first;
-      expr_t d1 = it.second;
-      code_file << FNUMEXPR_{ExpressionType::FirstExoDerivative, eq, var, lag};
-      if (prev_var != var || prev_lag != lag)
-        {
-          prev_var = var;
-          prev_lag = lag;
-          count_col_exo++;
-        }
-      d1->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-      code_file << FSTPG3_{eq, var, lag, count_col_exo-1};
-    }
-  // 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_{};
+                            simulation_type,
+                            0,
+                            symbol_table.endo_nbr(),
+                            endo_idx_block2orig,
+                            eq_idx_block2orig,
+                            false,
+                            symbol_table.endo_nbr(),
+                            max_endo_lag,
+                            max_endo_lead,
+                            u_count_int,
+                            jacobian_ncols_endo,
+                            symbol_table.exo_det_nbr(),
+                            jacobian_ncols_exo_det,
+                            symbol_table.exo_nbr(),
+                            jacobian_ncols_exo,
+                            0,
+                            0,
+                            exo_det,
+                            exo,
+                            {}};
+
+  writeBytecodeHelper<true>(code_file);
 }
 
 void
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index bef4a6bdd95792a63fb274ceef0c5d57a37ae5f7..4769d4923947f3fa903442d458e19487d550b44b 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1183,34 +1183,6 @@ ModelTree::testNestedParenthesis(const string &str) const
   return false;
 }
 
-void
-ModelTree::writeBytecodeTemporaryTerms(BytecodeWriter &code_file, ExprNodeBytecodeOutputType output_type, temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const
-{
-  // To store the functions that have already been written in the form TEF* = ext_fun();
-  for (auto [tt, idx] : temporary_terms_idxs)
-    {
-      if (dynamic_cast<AbstractExternalFunctionNode *>(tt))
-        tt->writeBytecodeExternalFunctionOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
-
-      code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, idx};
-      tt->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
-      switch (output_type)
-        {
-        case ExprNodeBytecodeOutputType::dynamicModel:
-          code_file << FSTPT_{idx};
-          break;
-        case ExprNodeBytecodeOutputType::staticModel:
-          code_file << FSTPST_{idx};
-          break;
-        case ExprNodeBytecodeOutputType::dynamicSteadyStateOperator:
-        case ExprNodeBytecodeOutputType::dynamicAssignmentLHS:
-        case ExprNodeBytecodeOutputType::staticAssignmentLHS:
-          cerr << "ModelTree::writeBytecodeTemporaryTerms: impossible case" << endl;
-          exit(EXIT_FAILURE);
-        }
-    }
-}
-
 void
 ModelTree::writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, deriv_node_temp_terms_t &tef_terms) const
 {
@@ -1255,79 +1227,39 @@ ModelTree::writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, d
   output << "]";
 }
 
-void
-ModelTree::writeBytecodeModelEquations(BytecodeWriter &code_file, ExprNodeBytecodeOutputType output_type, const temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
-{
-  for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
-    {
-      BinaryOpNode *eq_node = equations[eq];
-      expr_t lhs = eq_node->arg1, rhs = eq_node->arg2;
-      code_file << FNUMEXPR_{ExpressionType::ModelEquation, eq};
-      // Test if the right hand side of the equation is empty.
-      double vrhs = 1.0;
-      try
-        {
-          vrhs = rhs->eval({});
-        }
-      catch (ExprNode::EvalException &e)
-        {
-        }
-
-      if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
-        {
-          lhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
-          rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
-
-          code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{eq};
-        }
-      else // The right hand side of the equation is empty ==> residual=lhs;
-        {
-          lhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
-          code_file << FSTPR_{eq};
-        }
-    }
-}
-
-void
-ModelTree::writeBytecodeBinFile(const string &filename, int &u_count_int, bool &file_open,
-                                bool is_two_boundaries) const
+int
+ModelTree::writeBytecodeBinFile(const string &filename, bool is_two_boundaries) const
 {
-  int j;
-  std::ofstream SaveCode;
-  if (file_open)
-    SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate);
-  else
-    SaveCode.open(filename, ios::out | ios::binary);
+  ofstream SaveCode { 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;
-  for (const auto & [indices, d1] : derivatives[1])
-    {
-      int deriv_id = indices[1];
-      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
-        {
-          int eq = indices[0];
-          int var { getTypeSpecificIDByDerivID(deriv_id) };
-          int lag = getLagByDerivID(deriv_id);
-          SaveCode.write(reinterpret_cast<char *>(&eq), sizeof(eq));
-          int varr = var + lag * symbol_table.endo_nbr();
-          SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
-          SaveCode.write(reinterpret_cast<char *>(&lag), sizeof(lag));
-          int u = u_count_int + symbol_table.endo_nbr();
-          SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
-          u_count_int++;
-        }
-    }
+  int u_count {0};
+  for (const auto &[indices, d1] : derivatives[1])
+    if (int deriv_id {indices[1]};
+        getTypeByDerivID(deriv_id) == SymbolType::endogenous)
+      {
+        int eq {indices[0]};
+        SaveCode.write(reinterpret_cast<char *>(&eq), sizeof(eq));
+        int tsid {getTypeSpecificIDByDerivID(deriv_id)};
+        int lag {getLagByDerivID(deriv_id)};
+        int varr {tsid + lag * symbol_table.endo_nbr()};
+        SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
+        SaveCode.write(reinterpret_cast<char *>(&lag), sizeof(lag));
+        int u {u_count + symbol_table.endo_nbr()};
+        SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
+        u_count++;
+      }
   if (is_two_boundaries)
-    u_count_int += symbol_table.endo_nbr();
-  for (j = 0; j < symbol_table.endo_nbr(); j++)
+    u_count += symbol_table.endo_nbr();
+  for (int j {0}; j < symbol_table.endo_nbr(); j++)
     SaveCode.write(reinterpret_cast<char *>(&j), sizeof(j));
-  for (j = 0; j < symbol_table.endo_nbr(); j++)
+  for (int j {0}; j < symbol_table.endo_nbr(); j++)
     SaveCode.write(reinterpret_cast<char *>(&j), sizeof(j));
   SaveCode.close();
+  return u_count;
 }
 
 void
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index cdcdbd0183abcdbbb659812891285698391cf7a2..ad67aa1a7564e834f3e40e448af1b3d07518558b 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -242,9 +242,15 @@ protected:
   void writeTemporaryTerms(const temporary_terms_t &tt, temporary_terms_t &temp_term_union, const temporary_terms_idxs_t &tt_idxs, ostream &output, 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;
   //! Writes temporary terms in bytecode
-  void writeBytecodeTemporaryTerms(BytecodeWriter &code_file, ExprNodeBytecodeOutputType output_type, temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
-  //! Adds information for (non-block) bytecode simulation in a separate .bin file
-  void writeBytecodeBinFile(const string &filename, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
+  template<ExprNodeBytecodeOutputType output_type>
+  void writeBytecodeTemporaryTerms(const temporary_terms_t &tt,
+                                   temporary_terms_t &temporary_terms_union,
+                                   BytecodeWriter &code_file,
+                                   deriv_node_temp_terms_t &tef_terms) const;
+  /* Adds information for (non-block) bytecode simulation in a separate .bin
+     file.
+     Returns the number of first derivatives w.r.t. endogenous variables */
+  int writeBytecodeBinFile(const string &filename, bool is_two_boundaries) 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
@@ -269,6 +275,10 @@ protected:
   tuple<ostringstream, ostringstream, ostringstream, ostringstream,
         ostringstream, ostringstream, ostringstream> writeParamsDerivativesFileHelper() const;
 
+  // Helper for writing bytecode (without block decomposition)
+  template<bool dynamic>
+  void writeBytecodeHelper(BytecodeWriter &code_file) const;
+
   /* Helper for writing JSON output for residuals and derivatives.
      Returns mlv and derivatives output at each derivation order. */
   template<bool dynamic>
@@ -288,8 +298,10 @@ protected:
   /* Writes JSON model local variables.
      Optionally put the external function variable calls into TEF terms */
   void writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, deriv_node_temp_terms_t &tef_terms) const;
+
   //! Writes model equations in bytecode
-  void writeBytecodeModelEquations(BytecodeWriter &code_file, ExprNodeBytecodeOutputType output_type, const temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
+  template<ExprNodeBytecodeOutputType output_type>
+  void writeBytecodeModelEquations(BytecodeWriter &code_file, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms) const;
 
   //! Writes LaTeX model file
   void writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const;
@@ -1004,6 +1016,203 @@ ModelTree::writeParamsDerivativesFileHelper() const
     move(rpp_output), move(gpp_output), move(hp_output), move(g3p_output) };
 }
 
+template<ExprNodeBytecodeOutputType output_type>
+void
+ModelTree::writeBytecodeTemporaryTerms(const temporary_terms_t &tt,
+                                       temporary_terms_t &temporary_terms_union,
+                                       BytecodeWriter &code_file,
+                                       deriv_node_temp_terms_t &tef_terms) const
+{
+  // To store the functions that have already been written in the form TEF* = ext_fun();
+  for (auto it : tt)
+    {
+      if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+        it->writeBytecodeExternalFunctionOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
+
+      int idx {temporary_terms_idxs.at(it)};
+      code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, idx};
+      it->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
+      switch (output_type)
+        {
+        case ExprNodeBytecodeOutputType::dynamicModel:
+          code_file << FSTPT_{idx};
+          break;
+        case ExprNodeBytecodeOutputType::staticModel:
+          code_file << FSTPST_{idx};
+          break;
+        case ExprNodeBytecodeOutputType::dynamicSteadyStateOperator:
+        case ExprNodeBytecodeOutputType::dynamicAssignmentLHS:
+        case ExprNodeBytecodeOutputType::staticAssignmentLHS:
+          cerr << "ModelTree::writeBytecodeTemporaryTerms: impossible case" << endl;
+          exit(EXIT_FAILURE);
+
+          temporary_terms_union.insert(it);
+        }
+    }
+}
+
+template<ExprNodeBytecodeOutputType output_type>
+void
+ModelTree::writeBytecodeModelEquations(BytecodeWriter &code_file, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms) const
+{
+  for (int eq {0}; eq < static_cast<int>(equations.size()); eq++)
+    {
+      BinaryOpNode *eq_node {equations[eq]};
+      expr_t lhs {eq_node->arg1}, rhs {eq_node->arg2};
+      code_file << FNUMEXPR_{ExpressionType::ModelEquation, eq};
+      // Test if the right hand side of the equation is empty.
+      double vrhs {1.0};
+      try
+        {
+          vrhs = rhs->eval({});
+        }
+      catch (ExprNode::EvalException &e)
+        {
+        }
+
+      if (vrhs != 0) // The right hand side of the equation is not empty ⇒ residual=lhs-rhs
+        {
+          lhs->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+          rhs->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+
+          code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{eq};
+        }
+      else // The right hand side of the equation is empty ⇒ residual=lhs
+        {
+          lhs->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+          code_file << FSTPR_{eq};
+        }
+    }
+}
+
+template<bool dynamic>
+void
+ModelTree::writeBytecodeHelper(BytecodeWriter &code_file) const
+{
+  constexpr ExprNodeBytecodeOutputType output_type { dynamic ? ExprNodeBytecodeOutputType::dynamicModel : ExprNodeBytecodeOutputType::staticModel };
+
+  temporary_terms_t temporary_terms_union;
+  deriv_node_temp_terms_t tef_terms;
+
+  writeBytecodeTemporaryTerms<output_type>(temporary_terms_derivatives[0], temporary_terms_union, code_file, tef_terms);
+  writeBytecodeModelEquations<output_type>(code_file, temporary_terms_union, tef_terms);
+
+  code_file << FENDEQU_{};
+
+  // Temporary terms for the Jacobian
+  writeBytecodeTemporaryTerms<output_type>(temporary_terms_derivatives[1], temporary_terms_union, code_file, tef_terms);
+
+  // Get the current code_file position and jump if “evaluate” mode
+  int pos_jmpifeval {code_file.getInstructionCounter()};
+  code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being
+
+  // The Jacobian in “simulate” mode
+  vector<vector<tuple<int, int, int>>> my_derivatives(symbol_table.endo_nbr());;
+  int count_u {symbol_table.endo_nbr()};
+  for (const auto &[indices, d1] : derivatives[1])
+    {
+      auto [eq, deriv_id] {vectorToTuple<2>(indices)};
+      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
+        {
+          int tsid {getTypeSpecificIDByDerivID(deriv_id)};
+          int lag {getLagByDerivID(deriv_id)};
+          if constexpr(dynamic)
+            code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, tsid, lag};
+          else
+            code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, tsid};
+          if (!my_derivatives[eq].size())
+            my_derivatives[eq].clear();
+          my_derivatives[eq].emplace_back(tsid, lag, count_u);
+          d1->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
+          if constexpr(dynamic)
+            code_file << FSTPU_{count_u};
+          else
+            code_file << FSTPSU_{count_u};
+          count_u++;
+        }
+    }
+  for (int i {0}; i < symbol_table.endo_nbr(); i++)
+    {
+      code_file << FLDR_{i};
+      if (my_derivatives[i].size())
+        {
+          for (bool first_term {true};
+               const auto &[tsid, lag, uidx] : my_derivatives[i])
+            {
+              if constexpr(dynamic)
+                code_file << FLDU_{uidx} << FLDV_{SymbolType::endogenous, tsid, lag};
+              else
+                code_file << FLDSU_{uidx} << FLDSV_{SymbolType::endogenous, tsid};
+              code_file << FBINARY_{BinaryOpcode::times};
+              if (!exchange(first_term, false))
+                code_file << FBINARY_{BinaryOpcode::plus};
+            }
+          code_file << FBINARY_{BinaryOpcode::minus};
+        }
+      if constexpr(dynamic)
+        code_file << FSTPU_{i};
+      else
+        code_file << FSTPSU_{i};
+    }
+
+  // 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});
+
+  // The Jacobian in “evaluate” mode
+  for (const auto &[indices, d1] : derivatives[1])
+    {
+      auto [eq, deriv_id] {vectorToTuple<2>(indices)};
+      int tsid {getTypeSpecificIDByDerivID(deriv_id)};
+      int lag {getLagByDerivID(deriv_id)};
+      SymbolType type {getTypeByDerivID(deriv_id)};
+
+      if constexpr(dynamic)
+        {
+          ExpressionType expr_type;
+          switch (type)
+            {
+            case SymbolType::endogenous:
+              expr_type = ExpressionType::FirstEndoDerivative;
+              break;
+            case SymbolType::exogenous:
+              expr_type = ExpressionType::FirstExoDerivative;
+              break;
+            case SymbolType::exogenousDet:
+              expr_type = ExpressionType::FirstExodetDerivative;
+              break;
+            default:
+              assert(false);
+              break;
+            }
+          code_file << FNUMEXPR_{expr_type, eq, tsid, lag};
+        }
+      else
+        {
+          assert(type == SymbolType::endogenous);
+          code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, tsid};
+        }
+
+      d1->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs, tef_terms);
+      if constexpr(dynamic)
+         {
+           // Bytecode MEX uses a separate matrix for exogenous and exodet Jacobians
+           int jacob_col { type == SymbolType::endogenous ? getJacobianCol(deriv_id) : tsid };
+           code_file << FSTPG3_{eq, tsid, lag, jacob_col};
+         }
+      else
+        code_file << FSTPG2_{eq, tsid};
+    }
+
+  // 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_{};
+}
+
 template<bool dynamic>
 pair<ostringstream, vector<ostringstream>>
 ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 0d7af888a143643425330991619750df142c2be8..ad32e00720783e0e4724c22dd14e78a257d33074 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -261,118 +261,28 @@ StaticModel::writeStaticPerBlockCFiles(const string &basename) const
 void
 StaticModel::writeStaticBytecode(const string &basename) const
 {
-  ostringstream tmp_output;
-  bool file_open = false;
-
-  BytecodeWriter code_file{basename + "/model/bytecode/static.cod"};
-  int count_u;
-  int u_count_int = 0;
-
-  writeBytecodeBinFile(basename + "/model/bytecode/static.bin", u_count_int, file_open, false);
-  file_open = true;
-
-  // Compute the union of temporary terms from residuals and 1st derivatives
-  temporary_terms_t temporary_terms = temporary_terms_derivatives[0];
-  temporary_terms.insert(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end());
-
-  //Temporary variables declaration
-  code_file << FDIMST_{static_cast<int>(temporary_terms.size())}
-    << FBEGINBLOCK_{symbol_table.endo_nbr(),
-      BlockSimulationType::solveForwardComplete,
-      0,
-      symbol_table.endo_nbr(),
-      endo_idx_block2orig,
-      eq_idx_block2orig,
-      false,
-      symbol_table.endo_nbr(),
-      0,
-      0,
-      u_count_int,
-      symbol_table.endo_nbr()};
-
-  temporary_terms_t temporary_terms_union;
-  deriv_node_temp_terms_t tef_terms;
-
-  writeBytecodeTemporaryTerms(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-
-  writeBytecodeModelEquations(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-
-  code_file << FENDEQU_{};
-
-  // 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
-
-  vector<vector<pair<int, int>>> my_derivatives(symbol_table.endo_nbr());
-  count_u = symbol_table.endo_nbr();
-  for (const auto & [indices, d1] : derivatives[1])
-    {
-      int deriv_id = indices[1];
-      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
-        {
-          int eq = indices[0];
-          int var { getTypeSpecificIDByDerivID(deriv_id) };
-          code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, var};
-          if (!my_derivatives[eq].size())
-            my_derivatives[eq].clear();
-          my_derivatives[eq].emplace_back(var, count_u);
-
-          d1->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-
-          code_file << FSTPSU_{count_u};
-          count_u++;
-        }
-    }
-  for (int i = 0; i < symbol_table.endo_nbr(); i++)
-    {
-      code_file << FLDR_{i};
-      if (my_derivatives[i].size())
-        {
-          for (bool printed_something{false};
-               const auto &it : my_derivatives[i])
-            {
-              code_file << FLDSU_{it.second}
-                << FLDSV_{SymbolType::endogenous, it.first}
-                << FBINARY_{BinaryOpcode::times};
-              if (exchange(printed_something, true))
-                code_file << FBINARY_{BinaryOpcode::plus};
-            }
-          code_file << FBINARY_{BinaryOpcode::minus};
-        }
-      code_file << FSTPSU_{i};
-    }
-
-  // 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});
-
-  temporary_terms_t tt2, tt3;
-
-  // The Jacobian if we have to solve the block determinsitic bloc
-  for (const auto & [indices, d1] : derivatives[1])
-    {
-      int deriv_id = indices[1];
-      if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
-        {
-          int eq = indices[0];
-          int var { getTypeSpecificIDByDerivID(deriv_id) };
-          code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eq, var};
-          if (!my_derivatives[eq].size())
-            my_derivatives[eq].clear();
-          my_derivatives[eq].emplace_back(var, count_u);
-
-          d1->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, temporary_terms_idxs, tef_terms);
-          code_file << FSTPG2_{eq, var};
-        }
-    }
-
-  // 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_{};
+  // First write the .bin file
+  int u_count_int { writeBytecodeBinFile(basename + "/model/bytecode/static.bin", false) };
+
+  BytecodeWriter code_file {basename + "/model/bytecode/static.cod"};
+
+  // Declare temporary terms and the (single) block
+  code_file << FDIMST_{static_cast<int>(temporary_terms_derivatives[0].size()
+                                        + temporary_terms_derivatives[1].size())}
+            << FBEGINBLOCK_{symbol_table.endo_nbr(),
+                            BlockSimulationType::solveForwardComplete,
+                            0,
+                            symbol_table.endo_nbr(),
+                            endo_idx_block2orig,
+                            eq_idx_block2orig,
+                            false,
+                            symbol_table.endo_nbr(),
+                            0,
+                            0,
+                            u_count_int,
+                            symbol_table.endo_nbr()};
+
+  writeBytecodeHelper<false>(code_file);
 }
 
 void