diff --git a/src/Bytecode.cc b/src/Bytecode.cc
new file mode 100644
index 0000000000000000000000000000000000000000..510141697ec8d02eadf7b06b7cd6306cd3d5ff6f
--- /dev/null
+++ b/src/Bytecode.cc
@@ -0,0 +1,105 @@
+/*
+ * Copyright © 2022 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <iostream>
+#include <ios>
+#include <cstdlib>
+
+#include "Bytecode.hh"
+
+BytecodeWriter::BytecodeWriter(const string &filename)
+{
+  open(filename, ios::out | ios::binary);
+  if (!is_open())
+    {
+      cerr << R"(Error : Can't open file ")" << filename << R"(" for writing)" << endl;
+      exit(EXIT_FAILURE);
+    }
+}
+
+template<>
+BytecodeWriter &
+operator<<(BytecodeWriter &code_file, const FCALL_ &instr)
+{
+  code_file.instructions_positions.push_back(code_file.tellp());
+
+  code_file.write(reinterpret_cast<const char *>(&instr.op_code), sizeof instr.op_code);
+  code_file.write(reinterpret_cast<const char *>(&instr.nb_output_arguments), sizeof instr.nb_output_arguments);
+  code_file.write(reinterpret_cast<const char *>(&instr.nb_input_arguments), sizeof instr.nb_input_arguments);
+  code_file.write(reinterpret_cast<const char *>(&instr.indx), sizeof instr.indx);
+  code_file.write(reinterpret_cast<const char *>(&instr.add_input_arguments), sizeof instr.add_input_arguments);
+  code_file.write(reinterpret_cast<const char *>(&instr.row), sizeof instr.row);
+  code_file.write(reinterpret_cast<const char *>(&instr.col), sizeof instr.col);
+  code_file.write(reinterpret_cast<const char *>(&instr.function_type), sizeof instr.function_type);
+
+  int size = static_cast<int>(instr.func_name.size());
+  code_file.write(reinterpret_cast<char *>(&size), sizeof size);
+  const char *name = instr.func_name.c_str();
+  code_file.write(name, size);
+
+  size = static_cast<int>(instr.arg_func_name.size());
+  code_file.write(reinterpret_cast<char *>(&size), sizeof size);
+  name = instr.arg_func_name.c_str();
+  code_file.write(name, size);
+
+  return code_file;
+}
+
+template<>
+BytecodeWriter &
+operator<<(BytecodeWriter &code_file, const FBEGINBLOCK_ &instr)
+{
+  code_file.instructions_positions.push_back(code_file.tellp());
+
+  code_file.write(reinterpret_cast<const char *>(&instr.op_code), sizeof instr.op_code);
+  code_file.write(reinterpret_cast<const char *>(&instr.size), sizeof instr.size);
+  code_file.write(reinterpret_cast<const char *>(&instr.type), sizeof instr.type);
+  for (int i = 0; i < instr.size; i++)
+    {
+      code_file.write(reinterpret_cast<const char *>(&instr.variable[i]), sizeof instr.variable[i]);
+      code_file.write(reinterpret_cast<const char *>(&instr.equation[i]), sizeof instr.equation[i]);
+    }
+  if (instr.type == BlockSimulationType::solveTwoBoundariesSimple
+      || instr.type == BlockSimulationType::solveTwoBoundariesComplete
+      || instr.type == BlockSimulationType::solveBackwardComplete
+      || instr.type == BlockSimulationType::solveForwardComplete)
+    {
+      code_file.write(reinterpret_cast<const char *>(&instr.is_linear), sizeof instr.is_linear);
+      code_file.write(reinterpret_cast<const char *>(&instr.endo_nbr), sizeof instr.endo_nbr);
+      code_file.write(reinterpret_cast<const char *>(&instr.Max_Lag), sizeof instr.Max_Lag);
+      code_file.write(reinterpret_cast<const char *>(&instr.Max_Lead), sizeof instr.Max_Lead);
+      code_file.write(reinterpret_cast<const char *>(&instr.u_count_int), sizeof instr.u_count_int);
+    }
+  code_file.write(reinterpret_cast<const char *>(&instr.nb_col_jacob), sizeof instr.nb_col_jacob);
+  code_file.write(reinterpret_cast<const char *>(&instr.det_exo_size), sizeof instr.det_exo_size);
+  code_file.write(reinterpret_cast<const char *>(&instr.nb_col_det_exo_jacob), sizeof instr.nb_col_det_exo_jacob);
+  code_file.write(reinterpret_cast<const char *>(&instr.exo_size), sizeof instr.exo_size);
+  code_file.write(reinterpret_cast<const char *>(&instr.nb_col_exo_jacob), sizeof instr.nb_col_exo_jacob);
+  code_file.write(reinterpret_cast<const char *>(&instr.other_endo_size), sizeof instr.other_endo_size);
+  code_file.write(reinterpret_cast<const char *>(&instr.nb_col_other_endo_jacob), sizeof instr.nb_col_other_endo_jacob);
+
+  for (int i{0}; i < instr.det_exo_size; i++)
+    code_file.write(reinterpret_cast<const char *>(&instr.det_exogenous[i]), sizeof instr.det_exogenous[i]);
+  for (int i{0}; i < instr.exo_size; i++)
+    code_file.write(reinterpret_cast<const char *>(&instr.exogenous[i]), sizeof instr.exogenous[i]);
+  for (int i{0}; i < instr.other_endo_size; i++)
+    code_file.write(reinterpret_cast<const char *>(&instr.other_endogenous[i]), sizeof instr.other_endogenous[i]);
+
+  return code_file;
+}
diff --git a/src/Bytecode.hh b/src/Bytecode.hh
index ea4e56c3478e800ed728336f59547e9b1100be1d..6d179a5984147ef4339530ac184462d73f93c9dc 100644
--- a/src/Bytecode.hh
+++ b/src/Bytecode.hh
@@ -23,11 +23,13 @@
 #include <fstream>
 #include <vector>
 #include <utility>
+#include <ios>
+
+#include "CommonEnums.hh"
 
 #ifdef BYTECODE_MEX
 # include <dynmex.h>
 # include <cstring>
-# include "CommonEnums.hh"
 #endif
 
 using namespace std;
@@ -103,20 +105,18 @@ struct Block_contain_type
   int Equation, Variable, Own_Derivative;
 };
 
+class BytecodeWriter;
+
 class BytecodeInstruction
 {
+  template<typename B>
+  friend BytecodeWriter &operator<<(BytecodeWriter &code_file, const B &instr);
 protected:
   Tags op_code;
 public:
   explicit BytecodeInstruction(Tags op_code_arg) : op_code{op_code_arg}
   {
   };
-  void
-  write(ostream &CompileCode, unsigned int &instruction_number)
-  {
-    CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
-    instruction_number++;
-  };
 };
 
 template<typename T1>
@@ -129,12 +129,6 @@ public:
                                                       arg1{arg_arg1}
   {
   };
-  void
-  write(ostream &CompileCode, unsigned int &instruction_number)
-  {
-    CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
-    instruction_number++;
-  };
 };
 
 template<typename T1, typename T2>
@@ -148,12 +142,6 @@ public:
     BytecodeInstruction{op_code_arg}, arg1{arg_arg1}, arg2{arg_arg2}
   {
   };
-  void
-  write(ostream &CompileCode, unsigned int &instruction_number)
-  {
-    CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
-    instruction_number++;
-  };
 };
 
 template<typename T1, typename T2, typename T3>
@@ -168,12 +156,6 @@ public:
     BytecodeInstruction{op_code_arg}, arg1{arg_arg1}, arg2{arg_arg2}, arg3{arg_arg3}
   {
   };
-  void
-  write(ostream &CompileCode, unsigned int &instruction_number)
-  {
-    CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
-    instruction_number++;
-  };
 };
 
 template<typename T1, typename T2, typename T3, typename T4>
@@ -190,12 +172,6 @@ public:
     arg3{move(arg_arg3)}, arg4{arg_arg4}
   {
   };
-  void
-  write(ostream &CompileCode, unsigned int &instruction_number)
-  {
-    CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
-    instruction_number++;
-  };
 };
 
 class FLDZ_ : public BytecodeInstruction
@@ -750,6 +726,9 @@ public:
 
 class FCALL_ : public BytecodeInstruction
 {
+  template<typename B>
+  friend BytecodeWriter &operator<<(BytecodeWriter &code_file, const B &instr);
+private:
   int nb_output_arguments, nb_input_arguments, indx;
   string func_name;
   string arg_func_name;
@@ -838,27 +817,6 @@ public:
   {
     return function_type;
   }
-  void
-  write(ostream &CompileCode, unsigned int &instruction_number)
-  {
-    CompileCode.write(reinterpret_cast<char *>(&op_code), sizeof(op_code));
-    CompileCode.write(reinterpret_cast<char *>(&nb_output_arguments), sizeof(nb_output_arguments));
-    CompileCode.write(reinterpret_cast<char *>(&nb_input_arguments), sizeof(nb_input_arguments));
-    CompileCode.write(reinterpret_cast<char *>(&indx), sizeof(indx));
-    CompileCode.write(reinterpret_cast<char *>(&add_input_arguments), sizeof(add_input_arguments));
-    CompileCode.write(reinterpret_cast<char *>(&row), sizeof(row));
-    CompileCode.write(reinterpret_cast<char *>(&col), sizeof(col));
-    CompileCode.write(reinterpret_cast<char *>(&function_type), sizeof(function_type));
-    size_t size = func_name.size();
-    CompileCode.write(reinterpret_cast<char *>(&size), sizeof(int));
-    const char *name = func_name.c_str();
-    CompileCode.write(reinterpret_cast<const char *>(name), func_name.size());
-    size = arg_func_name.size();
-    CompileCode.write(reinterpret_cast<char *>(&size), sizeof(int));
-    name = arg_func_name.c_str();
-    CompileCode.write(reinterpret_cast<const char *>(name), arg_func_name.size());
-    instruction_number++;
-  };
 #ifdef BYTECODE_MEX
 
   char *
@@ -940,16 +898,12 @@ public:
   {
     return lag1;
   };
-  void
-  write(ostream &CompileCode, unsigned int &instruction_number)
-  {
-    CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
-    instruction_number++;
-  };
 };
 
 class FBEGINBLOCK_ : public BytecodeInstruction
 {
+  template<typename B>
+  friend BytecodeWriter &operator<<(BytecodeWriter &code_file, const B &instr);
 private:
   int size{0};
   BlockSimulationType type;
@@ -1106,44 +1060,6 @@ public:
   {
     return exogenous;
   }
-  void
-  write(ostream &CompileCode, unsigned int &instruction_number)
-  {
-    CompileCode.write(reinterpret_cast<char *>(&op_code), sizeof(op_code));
-    CompileCode.write(reinterpret_cast<char *>(&size), sizeof(size));
-    CompileCode.write(reinterpret_cast<char *>(&type), sizeof(type));
-    for (int i = 0; i < size; i++)
-      {
-        CompileCode.write(reinterpret_cast<char *>(&variable[i]), sizeof(variable[0]));
-        CompileCode.write(reinterpret_cast<char *>(&equation[i]), sizeof(equation[0]));
-      }
-    if (type == BlockSimulationType::solveTwoBoundariesSimple
-        || type == BlockSimulationType::solveTwoBoundariesComplete
-        || type == BlockSimulationType::solveBackwardComplete
-        || type == BlockSimulationType::solveForwardComplete)
-      {
-        CompileCode.write(reinterpret_cast<char *>(&is_linear), sizeof(is_linear));
-        CompileCode.write(reinterpret_cast<char *>(&endo_nbr), sizeof(endo_nbr));
-        CompileCode.write(reinterpret_cast<char *>(&Max_Lag), sizeof(Max_Lag));
-        CompileCode.write(reinterpret_cast<char *>(&Max_Lead), sizeof(Max_Lead));
-        CompileCode.write(reinterpret_cast<char *>(&u_count_int), sizeof(u_count_int));
-      }
-    CompileCode.write(reinterpret_cast<char *>(&nb_col_jacob), sizeof(nb_col_jacob));
-    CompileCode.write(reinterpret_cast<char *>(&det_exo_size), sizeof(det_exo_size));
-    CompileCode.write(reinterpret_cast<char *>(&nb_col_det_exo_jacob), sizeof(nb_col_det_exo_jacob));
-    CompileCode.write(reinterpret_cast<char *>(&exo_size), sizeof(exo_size));
-    CompileCode.write(reinterpret_cast<char *>(&nb_col_exo_jacob), sizeof(nb_col_exo_jacob));
-    CompileCode.write(reinterpret_cast<char *>(&other_endo_size), sizeof(other_endo_size));
-    CompileCode.write(reinterpret_cast<char *>(&nb_col_other_endo_jacob), sizeof(nb_col_other_endo_jacob));
-
-    for (int i{0}; i < det_exo_size; i++)
-      CompileCode.write(reinterpret_cast<char *>(&det_exogenous[i]), sizeof(det_exogenous[0]));
-    for (int i{0}; i < exo_size; i++)
-      CompileCode.write(reinterpret_cast<char *>(&exogenous[i]), sizeof(exogenous[0]));
-    for (int i{0}; i < other_endo_size; i++)
-      CompileCode.write(reinterpret_cast<char *>(&other_endogenous[i]), sizeof(other_endogenous[0]));
-    instruction_number++;
-  };
 #ifdef BYTECODE_MEX
 
   char *
@@ -1201,6 +1117,54 @@ public:
 #endif
 };
 
+// Superclass of std::ofstream for writing a sequence of bytecode instructions
+class BytecodeWriter : private ofstream
+{
+  template<typename B>
+  friend BytecodeWriter &operator<<(BytecodeWriter &code_file, const B &instr);
+private:
+  // Stores the positions of all instructions in the byte stream
+  vector<pos_type> instructions_positions;
+public:
+  BytecodeWriter(const string &filename);
+  // Returns the number of the next instruction to be written
+  int
+  getInstructionCounter() const
+  {
+    return static_cast<int>(instructions_positions.size());
+  }
+  /* Overwrites an existing instruction, given its number.
+     It is the responsibility of the caller to ensure that the new instruction
+     occupies exactly as many bytes as the former one. */
+  template<typename B>
+  void
+  overwriteInstruction(int instruction_number, const B &new_instruction)
+  {
+    seekp(instructions_positions.at(instruction_number));
+    *this << new_instruction;
+    instructions_positions.pop_back();
+    seekp(0, ios_base::end);
+  }
+};
+
+// Overloads of operator<< for writing bytecode instructions
+
+template<typename B>
+BytecodeWriter &
+operator<<(BytecodeWriter &code_file, const B &instr)
+{
+  code_file.instructions_positions.push_back(code_file.tellp());
+  code_file.write(reinterpret_cast<const char *>(&instr), sizeof(B));
+  return code_file;
+}
+
+template<>
+BytecodeWriter &operator<<(BytecodeWriter &code_file, const FCALL_ &instr);
+
+template<>
+BytecodeWriter &operator<<(BytecodeWriter &code_file, const FBEGINBLOCK_ &instr);
+
+
 #ifdef BYTECODE_MEX
 using tags_liste_t = vector<pair<Tags, void * >>;
 class CodeLoad
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 7225700120aa460b12ab2a9a57c6d2d8edd5d157..e165fe2064fb62a647d4e535e95978417311ca99 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -28,7 +28,6 @@
 
 #include "DynamicModel.hh"
 #include "ParsingDriver.hh"
-#include "Bytecode.hh"
 
 void
 DynamicModel::copyHelper(const DynamicModel &m)
@@ -169,29 +168,23 @@ DynamicModel::operator=(const DynamicModel &m)
 }
 
 void
-DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, 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
+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->compile(code_file, instruction_number, false, temporary_terms, temporary_terms_idxs, true, false, tef_terms);
+    it->second->writeBytecodeOutput(code_file, false, temporary_terms, temporary_terms_idxs, true, false, tef_terms);
   else
-    {
-      FLDZ_ fldz;
-      fldz.write(code_file, instruction_number);
-    }
+    code_file << FLDZ_{};
 }
 
 void
-DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, 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
+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->compile(code_file, instruction_number, false, temporary_terms, temporary_terms_idxs, true, false, tef_terms);
+    it->second->writeBytecodeOutput(code_file, false, temporary_terms, temporary_terms_idxs, true, false, tef_terms);
   else
-    {
-      FLDZ_ fldz;
-      fldz.write(code_file, instruction_number);
-    }
+    code_file << FLDZ_{};
 }
 
 void
@@ -776,18 +769,9 @@ void
 DynamicModel::writeDynamicBytecode(const string &basename) const
 {
   ostringstream tmp_output;
-  ofstream code_file;
-  unsigned int instruction_number = 0;
+  BytecodeWriter code_file{basename + "/model/bytecode/dynamic.cod"};
   bool file_open = false;
 
-  string main_name = basename + "/model/bytecode/dynamic.cod";
-  code_file.open(main_name, ios::out | ios::binary | ios::ate);
-  if (!code_file.is_open())
-    {
-      cerr << R"(Error : Can't open file ")" << main_name << R"(" for writing)" << endl;
-      exit(EXIT_FAILURE);
-    }
-
   int count_u;
   int u_count_int = 0;
   BlockSimulationType simulation_type;
@@ -802,8 +786,7 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
   file_open = true;
 
   //Temporary variables declaration
-  FDIMT_ fdimt{static_cast<int>(temporary_terms_idxs.size())};
-  fdimt.write(code_file, instruction_number);
+  code_file << FDIMT_{static_cast<int>(temporary_terms_idxs.size())};
 
   vector<int> exo, exo_det, other_endo;
 
@@ -863,44 +846,40 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
         }
     }
 
-  FBEGINBLOCK_ 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);
-  fbeginblock.write(code_file, instruction_number);
+  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;
 
-  compileTemporaryTerms(code_file, instruction_number, true, false, temporary_terms_union, temporary_terms_idxs, tef_terms);
+  writeBytecodeTemporaryTerms(code_file, true, false, temporary_terms_union, temporary_terms_idxs, tef_terms);
 
-  compileModelEquations(code_file, instruction_number, true, false, temporary_terms_union, temporary_terms_idxs, tef_terms);
+  writeBytecodeModelEquations(code_file, true, false, temporary_terms_union, temporary_terms_idxs, tef_terms);
 
-  FENDEQU_ fendequ;
-  fendequ.write(code_file, instruction_number);
+  code_file << FENDEQU_{};
 
-  // Get the current code_file position and jump if eval = true
-  streampos pos1 = code_file.tellp();
-  FJMPIFEVAL_ fjmp_if_eval(0);
-  fjmp_if_eval.write(code_file, instruction_number);
-  int prev_instruction_number = instruction_number;
+  // 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();
@@ -913,56 +892,39 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
           int symb = getSymbIDByDerivID(deriv_id);
           int var = symbol_table.getTypeSpecificID(symb);
           int lag = getLagByDerivID(deriv_id);
-          FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eq, var, lag);
-          fnumexpr.write(code_file, instruction_number);
+          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->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, true, false, tef_terms);
+          d1->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, true, false, tef_terms);
 
-          FSTPU_ fstpu(count_u);
-          fstpu.write(code_file, instruction_number);
+          code_file << FSTPU_{count_u};
           count_u++;
         }
     }
   for (int i = 0; i < symbol_table.endo_nbr(); i++)
     {
-      FLDR_ fldr(i);
-      fldr.write(code_file, instruction_number);
+      code_file << FLDR_{i};
       if (my_derivatives[i].size())
         {
           for (auto it = my_derivatives[i].begin(); it != my_derivatives[i].end(); ++it)
             {
-              FLDU_ fldu(get<2>(*it));
-              fldu.write(code_file, instruction_number);
-              FLDV_ fldv{SymbolType::endogenous, get<0>(*it), get<1>(*it)};
-              fldv.write(code_file, instruction_number);
-              FBINARY_ fbinary{BinaryOpcode::times};
-              fbinary.write(code_file, instruction_number);
+              code_file << FLDU_{get<2>(*it)}
+                << FLDV_{SymbolType::endogenous, get<0>(*it), get<1>(*it)}
+                << FBINARY_{BinaryOpcode::times};
               if (it != my_derivatives[i].begin())
-                {
-                  FBINARY_ fbinary{BinaryOpcode::plus};
-                  fbinary.write(code_file, instruction_number);
-                }
+                code_file << FBINARY_{BinaryOpcode::plus};
             }
-          FBINARY_ fbinary{BinaryOpcode::minus};
-          fbinary.write(code_file, instruction_number);
+          code_file << FBINARY_{BinaryOpcode::minus};
         }
-      FSTPU_ fstpu(i);
-      fstpu.write(code_file, instruction_number);
-    }
-
-  // Get the current code_file position and jump = true
-  streampos pos2 = code_file.tellp();
-  FJMP_ fjmp(0);
-  fjmp.write(code_file, instruction_number);
-  // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
-  streampos pos3 = code_file.tellp();
-  code_file.seekp(pos1);
-  FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
-  fjmp_if_eval1.write(code_file, instruction_number);
-  code_file.seekp(pos3);
-  prev_instruction_number = instruction_number;
+      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;
@@ -972,17 +934,15 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
     {
       auto [lag, var, eq] = it.first;
       expr_t d1 = it.second;
-      FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eq, var, lag);
-      fnumexpr.write(code_file, instruction_number);
+      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->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, true, false, tef_terms);
-      FSTPG3_ fstpg3(eq, var, lag, count_col_endo-1);
-      fstpg3.write(code_file, instruction_number);
+      d1->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, true, false, tef_terms);
+      code_file << FSTPG3_{eq, var, lag, count_col_endo-1};
     }
   prev_var = -1;
   prev_lag = -999999999;
@@ -991,30 +951,21 @@ DynamicModel::writeDynamicBytecode(const string &basename) const
     {
       auto [lag, ignore, var, eq] = it.first;
       expr_t d1 = it.second;
-      FNUMEXPR_ fnumexpr(ExpressionType::FirstExoDerivative, eq, var, lag);
-      fnumexpr.write(code_file, instruction_number);
+      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->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, true, false, tef_terms);
-      FSTPG3_ fstpg3(eq, var, lag, count_col_exo-1);
-      fstpg3.write(code_file, instruction_number);
-    }
-  // Set codefile position to previous JMP_ and set the number of instructions to jump
-  pos1 = code_file.tellp();
-  code_file.seekp(pos2);
-  FJMP_ fjmp1(instruction_number - prev_instruction_number);
-  fjmp1.write(code_file, instruction_number);
-  code_file.seekp(pos1);
-
-  FENDBLOCK_ fendblock;
-  fendblock.write(code_file, instruction_number);
-  FEND_ fend;
-  fend.write(code_file, instruction_number);
-  code_file.close();
+      d1->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, true, false, 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_{};
 }
 
 void
@@ -1034,34 +985,23 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
   int i, v;
   string tmp_s;
   ostringstream tmp_output;
-  ofstream code_file;
-  unsigned int instruction_number = 0;
   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;
-  string main_name = basename + "/model/bytecode/dynamic.cod";
-  code_file.open(main_name, ios::out | ios::binary | ios::ate);
-  if (!code_file.is_open())
-    {
-      cerr << R"(Error : Can't open file ")" << main_name << R"(" for writing)" << endl;
-      exit(EXIT_FAILURE);
-    }
+  BytecodeWriter code_file{basename + "/model/bytecode/dynamic.cod"};
+
   //Temporary variables declaration
 
-  FDIMT_ fdimt{static_cast<int>(blocks_temporary_terms_idxs.size())};
-  fdimt.write(code_file, instruction_number);
+  code_file << FDIMT_{static_cast<int>(blocks_temporary_terms_idxs.size())};
 
   for (int block = 0; block < static_cast<int>(blocks.size()); block++)
     {
       feedback_variables.clear();
       if (block > 0)
-        {
-          FENDBLOCK_ fendblock;
-          fendblock.write(code_file, instruction_number);
-        }
+        code_file << FENDBLOCK_{};
       int count_u;
       int u_count_int = 0;
       BlockSimulationType simulation_type = blocks[block].simulation_type;
@@ -1081,28 +1021,27 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
           file_open = true;
         }
 
-      FBEGINBLOCK_ 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,
-                               blocks_jacob_cols_endo[block].size(),
-                               blocks_exo_det[block].size(),
-                               blocks_jacob_cols_exo_det[block].size(),
-                               blocks_exo[block].size(),
-                               blocks_jacob_cols_exo[block].size(),
-                               blocks_other_endo[block].size(),
-                               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()));
-      fbeginblock.write(code_file, instruction_number);
+      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;
 
@@ -1114,13 +1053,11 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
                            for (auto it : blocks_temporary_terms[block][eq])
                              {
                                if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-                                 it->compileExternalFunctionOutput(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+                                 it->writeBytecodeExternalFunctionOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
 
-                               FNUMEXPR_ fnumexpr{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)};
-                               fnumexpr.write(code_file, instruction_number);
-                               it->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
-                               FSTPT_ fstpt{blocks_temporary_terms_idxs.at(it)};
-                               fstpt.write(code_file, instruction_number);
+                               code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)};
+                               it->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+                               code_file << FSTPT_{blocks_temporary_terms_idxs.at(it)};
                                temporary_terms_union.insert(it);
                              }
                          };
@@ -1139,25 +1076,22 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
             case BlockSimulationType::evaluateBackward:
             case BlockSimulationType::evaluateForward:
               equ_type = getBlockEquationType(block, i);
-              {
-                FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, getBlockEquationID(block, i));
-                fnumexpr.write(code_file, instruction_number);
-              }
+              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->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+                  rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+                  lhs->writeBytecodeOutput(code_file, true, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
                 }
               else if (equ_type == EquationType::evaluateRenormalized)
                 {
                   eq_node = getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+                  rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+                  lhs->writeBytecodeOutput(code_file, true, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
                 }
               break;
             case BlockSimulationType::solveBackwardComplete:
@@ -1173,22 +1107,17 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
               goto end;
             default:
             end:
-              FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, getBlockEquationID(block, i));
-              fnumexpr.write(code_file, instruction_number);
+              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
               eq_node = getBlockEquationExpr(block, i);
               lhs = eq_node->arg1;
               rhs = eq_node->arg2;
-              lhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
-              rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+              lhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+              rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
 
-              FBINARY_ fbinary{BinaryOpcode::minus};
-              fbinary.write(code_file, instruction_number);
-              FSTPR_ fstpr(i - block_recursive);
-              fstpr.write(code_file, instruction_number);
+              code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
             }
         }
-      FENDEQU_ fendequ;
-      fendequ.write(code_file, instruction_number);
+      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
@@ -1198,11 +1127,9 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
           && simulation_type != BlockSimulationType::evaluateForward)
         write_eq_tt(blocks[block].size);
 
-      // Get the current code_file position and jump if eval = true
-      streampos pos1 = code_file.tellp();
-      FJMPIFEVAL_ fjmp_if_eval(0);
-      fjmp_if_eval.write(code_file, instruction_number);
-      int prev_instruction_number = instruction_number;
+      // 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”) */
@@ -1213,15 +1140,9 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
             {
             case BlockSimulationType::solveBackwardSimple:
             case BlockSimulationType::solveForwardSimple:
-              {
-                FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0);
-                fnumexpr.write(code_file, instruction_number);
-              }
-              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-              {
-                FSTPG_ fstpg(0);
-                fstpg.write(code_file, instruction_number);
-              }
+              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:
@@ -1254,11 +1175,9 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
                       Uf[eqr].Ufl->u = count_u;
                       Uf[eqr].Ufl->var = varr;
                       Uf[eqr].Ufl->lag = lag;
-                      FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eqr, varr, lag);
-                      fnumexpr.write(code_file, instruction_number);
-                      compileChainRuleDerivative(code_file, instruction_number, block, eq, var, lag, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                      FSTPU_ fstpu(count_u);
-                      fstpu.write(code_file, instruction_number);
+                      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++;
                     }
                 }
@@ -1266,26 +1185,14 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
                 {
                   if (i >= block_recursive)
                     {
-                      FLDR_ fldr(i-block_recursive);
-                      fldr.write(code_file, instruction_number);
-
-                      FLDZ_ fldz;
-                      fldz.write(code_file, instruction_number);
+                      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)
-                        {
-                          FLDU_ fldu(Uf[v].Ufl->u);
-                          fldu.write(code_file, instruction_number);
-                          FLDV_ fldv{SymbolType::endogenous, Uf[v].Ufl->var, Uf[v].Ufl->lag};
-                          fldv.write(code_file, instruction_number);
-
-                          FBINARY_ fbinary{BinaryOpcode::times};
-                          fbinary.write(code_file, instruction_number);
-
-                          FCUML_ fcuml;
-                          fcuml.write(code_file, instruction_number);
-                        }
+                        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)
                         {
@@ -1293,11 +1200,8 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
                           free(Uf[v].Ufl);
                           Uf[v].Ufl = Uf[v].Ufl_First;
                         }
-                      FBINARY_ fbinary{BinaryOpcode::minus};
-                      fbinary.write(code_file, instruction_number);
-
-                      FSTPU_ fstpu(i - block_recursive);
-                      fstpu.write(code_file, instruction_number);
+                      code_file << FBINARY_{BinaryOpcode::minus}
+                        << FSTPU_{i - block_recursive};
                     }
                 }
               break;
@@ -1306,17 +1210,11 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
             }
         }
 
-      // Get the current code_file position and jump = true
-      streampos pos2 = code_file.tellp();
-      FJMP_ fjmp(0);
-      fjmp.write(code_file, instruction_number);
-      // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
-      streampos pos3 = code_file.tellp();
-      code_file.seekp(pos1);
-      FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
-      fjmp_if_eval1.write(code_file, instruction_number);
-      code_file.seekp(pos3);
-      prev_instruction_number = instruction_number;
+      // 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
@@ -1331,58 +1229,43 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const
           auto [eq, var, lag] = indices;
           int eqr = getBlockEquationID(block, eq);
           int varr = getBlockVariableID(block, var);
-          FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eqr, varr, lag);
-          fnumexpr.write(code_file, instruction_number);
-          compileDerivative(code_file, instruction_number, eqr, varr, lag, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-          FSTPG3_ fstpg3(eq, var, lag, blocks_jacob_cols_endo[block].at({ var, lag }));
-          fstpg3.write(code_file, instruction_number);
+          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
-          FNUMEXPR_ fnumexpr(ExpressionType::FirstExoDerivative, eqr, varr, lag);
-          fnumexpr.write(code_file, instruction_number);
-          d->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
-          FSTPG3_ fstpg3(eq, var, lag, blocks_jacob_cols_exo[block].at({ var, lag }));
-          fstpg3.write(code_file, instruction_number);
+          code_file << FNUMEXPR_{ExpressionType::FirstExoDerivative, eqr, varr, lag};
+          d->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, 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
-          FNUMEXPR_ fnumexpr(ExpressionType::FirstExodetDerivative, eqr, varr, lag);
-          fnumexpr.write(code_file, instruction_number);
-          d->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
-          FSTPG3_ fstpg3(eq, var, lag, blocks_jacob_cols_exo_det[block].at({ var, lag }));
-          fstpg3.write(code_file, instruction_number);
+          code_file << FNUMEXPR_{ExpressionType::FirstExodetDerivative, eqr, varr, lag};
+          d->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, 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
-          FNUMEXPR_ fnumexpr(ExpressionType::FirstOtherEndoDerivative, eqr, varr, lag);
-          fnumexpr.write(code_file, instruction_number);
-          d->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
-          FSTPG3_ fstpg3(eq, var, lag, blocks_jacob_cols_other_endo[block].at({ var, lag }));
-          fstpg3.write(code_file, instruction_number);
+          code_file << FNUMEXPR_{ExpressionType::FirstOtherEndoDerivative, eqr, varr, lag};
+          d->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+          code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_other_endo[block].at({ var, lag })};
         }
 
-      // Set codefile position to previous JMP_ and set the number of instructions to jump
-      pos1 = code_file.tellp();
-      code_file.seekp(pos2);
-      FJMP_ fjmp1(instruction_number - prev_instruction_number);
-      fjmp1.write(code_file, instruction_number);
-      code_file.seekp(pos1);
-    }
-  FENDBLOCK_ fendblock;
-  fendblock.write(code_file, instruction_number);
-  FEND_ fend;
-  fend.write(code_file, instruction_number);
-  code_file.close();
+      // 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_{};
 }
 
 void
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 191683152a4a49203d8603d14c74f28fd034a2e7..0d5d96aeda6fb1ad9d40c07cea537099a94f3f72 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -24,6 +24,7 @@
 #include <filesystem>
 
 #include "StaticModel.hh"
+#include "Bytecode.hh"
 
 using namespace std;
 
@@ -182,10 +183,10 @@ private:
                                      vector<vector<temporary_terms_t>> &blocks_temporary_terms,
                                      map<expr_t, tuple<int, int, int>> &reference_count) const override;
 
-  //! Write derivative code of an equation w.r. to a variable
-  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const 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 code of an equation w.r. to a variable
-  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
+  //! 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;
 
   //! Get the type corresponding to a derivation ID
   SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 7038847407275afe79d499b2227c2c6cebd12a74..8216988ed42385ec617132974406db01b92d9226 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -27,7 +27,6 @@
 #include "ExprNode.hh"
 #include "DataTree.hh"
 #include "ModFile.hh"
-#include "Bytecode.hh"
 
 ExprNode::ExprNode(DataTree &datatree_arg, int idx_arg) : datatree{datatree_arg}, idx{idx_arg}
 {
@@ -107,6 +106,28 @@ ExprNode::checkIfTemporaryTermThenWrite(ostream &output, ExprNodeOutputType outp
   return true;
 }
 
+bool
+ExprNode::checkIfTemporaryTermThenWriteBytecode(BytecodeWriter &code_file,
+                                                const temporary_terms_t &temporary_terms,
+                                                const temporary_terms_idxs_t &temporary_terms_idxs,
+                                                bool dynamic) const
+{
+  if (!temporary_terms.contains(const_cast<ExprNode *>(this)))
+    return false;
+
+  auto it2 = temporary_terms_idxs.find(const_cast<ExprNode *>(this));
+  // It is the responsibility of the caller to ensure that all temporary terms have their index
+  assert(it2 != temporary_terms_idxs.end());
+
+  if (dynamic)
+    code_file << FLDT_{it2->second};
+  else
+    code_file << FLDST_{it2->second};
+
+  return true;
+}
+
+
 pair<expr_t, int>
 ExprNode::getLagEquivalenceClass() const
 {
@@ -189,10 +210,10 @@ ExprNode::writeJsonExternalFunctionOutput(vector<string> &efout,
 }
 
 void
-ExprNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                        const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                        deriv_node_temp_terms_t &tef_terms) const
+ExprNode::writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                              const temporary_terms_t &temporary_terms,
+                                              const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                              deriv_node_temp_terms_t &tef_terms) const
 {
   // Nothing to do
 }
@@ -456,13 +477,12 @@ NumConstNode::eval(const eval_context_t &eval_context) const noexcept(false)
 }
 
 void
-NumConstNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                      const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                      const deriv_node_temp_terms_t &tef_terms) const
+NumConstNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                  const temporary_terms_t &temporary_terms,
+                                  const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                  const deriv_node_temp_terms_t &tef_terms) const
 {
-  FLDC_ fldc(datatree.num_constants.getDouble(id));
-  fldc.write(CompileCode, instruction_number);
+  code_file << FLDC_{datatree.num_constants.getDouble(id)};
 }
 
 void
@@ -1225,14 +1245,14 @@ VariableNode::eval(const eval_context_t &eval_context) const noexcept(false)
 }
 
 void
-VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                      const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                      const deriv_node_temp_terms_t &tef_terms) const
+VariableNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                  const temporary_terms_t &temporary_terms,
+                                  const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                  const deriv_node_temp_terms_t &tef_terms) const
 {
   auto type = get_type();
   if (type == SymbolType::modelLocalVariable || type == SymbolType::modFileLocalVariable)
-    datatree.getLocalVariable(symb_id)->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+    datatree.getLocalVariable(symb_id)->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
   else
     {
       int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
@@ -1243,29 +1263,17 @@ VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number,
           if (dynamic)
             {
               if (steady_dynamic) // steady state values in a dynamic model
-                {
-                  FLDVS_ fldvs{type, tsid};
-                  fldvs.write(CompileCode, instruction_number);
-                }
+                code_file << FLDVS_{type, tsid};
               else
                 {
                   if (type == SymbolType::parameter)
-                    {
-                      FLDV_ fldv{type, tsid};
-                      fldv.write(CompileCode, instruction_number);
-                    }
+                    code_file << FLDV_{type, tsid};
                   else
-                    {
-                      FLDV_ fldv{type, tsid, lag};
-                      fldv.write(CompileCode, instruction_number);
-                    }
+                    code_file << FLDV_{type, tsid, lag};
                 }
             }
           else
-            {
-              FLDSV_ fldsv{type, tsid};
-              fldsv.write(CompileCode, instruction_number);
-            }
+            code_file << FLDSV_{type, tsid};
         }
       else
         {
@@ -1279,22 +1287,13 @@ VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number,
               else
                 {
                   if (type == SymbolType::parameter)
-                    {
-                      FSTPV_ fstpv{type, tsid};
-                      fstpv.write(CompileCode, instruction_number);
-                    }
+                    code_file << FSTPV_{type, tsid};
                   else
-                    {
-                      FSTPV_ fstpv{type, tsid, lag};
-                      fstpv.write(CompileCode, instruction_number);
-                    }
+                    code_file << FSTPV_{type, tsid, lag};
                 }
             }
           else
-            {
-              FSTPSV_ fstpsv{type, tsid};
-              fstpsv.write(CompileCode, instruction_number);
-            }
+            code_file << FSTPSV_{type, tsid};
         }
     }
 }
@@ -2923,13 +2922,13 @@ UnaryOpNode::writeJsonExternalFunctionOutput(vector<string> &efout,
 }
 
 void
-UnaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                           deriv_node_temp_terms_t &tef_terms) const
+UnaryOpNode::writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                                 const temporary_terms_t &temporary_terms,
+                                                 const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                 deriv_node_temp_terms_t &tef_terms) const
 {
-  arg->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                     temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs, temporary_terms,
+                                           temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 }
 
 double
@@ -3011,33 +3010,20 @@ UnaryOpNode::eval(const eval_context_t &eval_context) const noexcept(false)
 }
 
 void
-UnaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                     const deriv_node_temp_terms_t &tef_terms) const
+UnaryOpNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                 const temporary_terms_t &temporary_terms,
+                                 const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                 const deriv_node_temp_terms_t &tef_terms) const
 {
-  if (auto this2 = const_cast<UnaryOpNode *>(this);
-      temporary_terms.contains(this2))
-    {
-      if (dynamic)
-        {
-          FLDT_ fldt(temporary_terms_idxs.at(this2));
-          fldt.write(CompileCode, instruction_number);
-        }
-      else
-        {
-          FLDST_ fldst(temporary_terms_idxs.at(this2));
-          fldst.write(CompileCode, instruction_number);
-        }
-      return;
-    }
+  if (checkIfTemporaryTermThenWriteBytecode(code_file, temporary_terms, temporary_terms_idxs, dynamic))
+    return;
+
   if (op_code == UnaryOpcode::steadyState)
-    arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, true, tef_terms);
+    arg->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, true, tef_terms);
   else
     {
-      arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-      FUNARY_ funary{op_code};
-      funary.write(CompileCode, instruction_number);
+      arg->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+      code_file << FUNARY_{op_code};
     }
 }
 
@@ -4271,36 +4257,19 @@ BinaryOpNode::eval(const eval_context_t &eval_context) const noexcept(false)
 }
 
 void
-BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                      const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                      const deriv_node_temp_terms_t &tef_terms) const
+BinaryOpNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                  const temporary_terms_t &temporary_terms,
+                                  const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                  const deriv_node_temp_terms_t &tef_terms) const
 {
-  // If current node is a temporary term
-  if (auto this2 = const_cast<BinaryOpNode *>(this);
-      temporary_terms.contains(this2))
-    {
-      if (dynamic)
-        {
-          FLDT_ fldt(temporary_terms_idxs.at(this2));
-          fldt.write(CompileCode, instruction_number);
-        }
-      else
-        {
-          FLDST_ fldst(temporary_terms_idxs.at(this2));
-          fldst.write(CompileCode, instruction_number);
-        }
-      return;
-    }
+  if (checkIfTemporaryTermThenWriteBytecode(code_file, temporary_terms, temporary_terms_idxs, dynamic))
+    return;
+
   if (op_code == BinaryOpcode::powerDeriv)
-    {
-      FLDC_ fldc(powerDerivOrder);
-      fldc.write(CompileCode, instruction_number);
-    }
-  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-  FBINARY_ fbinary{op_code};
-  fbinary.write(CompileCode, instruction_number);
+    code_file << FLDC_{static_cast<double>(powerDerivOrder)};
+  arg1->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg2->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  code_file << FBINARY_{op_code};
 }
 
 bool
@@ -4738,15 +4707,15 @@ BinaryOpNode::writeJsonExternalFunctionOutput(vector<string> &efout,
 }
 
 void
-BinaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                            bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                            const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                            deriv_node_temp_terms_t &tef_terms) const
+BinaryOpNode::writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                                  const temporary_terms_t &temporary_terms,
+                                                  const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                  deriv_node_temp_terms_t &tef_terms) const
 {
-  arg1->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                      temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-  arg2->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                      temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg1->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs, temporary_terms,
+                                            temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg2->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs, temporary_terms,
+                                            temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 }
 
 int
@@ -5971,32 +5940,18 @@ TrinaryOpNode::eval(const eval_context_t &eval_context) const noexcept(false)
 }
 
 void
-TrinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                       bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                       const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                       const deriv_node_temp_terms_t &tef_terms) const
+TrinaryOpNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                   const temporary_terms_t &temporary_terms,
+                                   const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                   const deriv_node_temp_terms_t &tef_terms) const
 {
-  // If current node is a temporary term
-  if (auto this2 = const_cast<TrinaryOpNode *>(this);
-      temporary_terms.contains(this2))
-    {
-      if (dynamic)
-        {
-          FLDT_ fldt(temporary_terms_idxs.at(this2));
-          fldt.write(CompileCode, instruction_number);
-        }
-      else
-        {
-          FLDST_ fldst(temporary_terms_idxs.at(this2));
-          fldst.write(CompileCode, instruction_number);
-        }
-      return;
-    }
-  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-  arg3->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-  FTRINARY_ ftrinary{op_code};
-  ftrinary.write(CompileCode, instruction_number);
+  if (checkIfTemporaryTermThenWriteBytecode(code_file, temporary_terms, temporary_terms_idxs, dynamic))
+    return;
+
+  arg1->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg2->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg3->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  code_file << FTRINARY_{op_code};
 }
 
 bool
@@ -6174,17 +6129,17 @@ TrinaryOpNode::writeJsonExternalFunctionOutput(vector<string> &efout,
 }
 
 void
-TrinaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                             bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                             const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                             deriv_node_temp_terms_t &tef_terms) const
+TrinaryOpNode::writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                                   const temporary_terms_t &temporary_terms,
+                                                   const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                   deriv_node_temp_terms_t &tef_terms) const
 {
-  arg1->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                      temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-  arg2->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                      temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-  arg3->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                      temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg1->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs, temporary_terms,
+                                            temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg2->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs, temporary_terms,
+                                            temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+  arg3->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs, temporary_terms,
+                                            temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 }
 
 void
@@ -6655,14 +6610,14 @@ AbstractExternalFunctionNode::getChainRuleDerivative(int deriv_id, const map<int
 }
 
 int
-AbstractExternalFunctionNode::compileExternalFunctionArguments(ostream &CompileCode, unsigned int &instruction_number,
-                                                               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                                               const deriv_node_temp_terms_t &tef_terms) const
+AbstractExternalFunctionNode::writeBytecodeExternalFunctionArguments(BytecodeWriter &code_file, bool lhs_rhs,
+                                                                     const temporary_terms_t &temporary_terms,
+                                                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                                     const deriv_node_temp_terms_t &tef_terms) const
 {
   for (auto argument : arguments)
-    argument->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                      temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+    argument->writeBytecodeOutput(code_file, lhs_rhs, temporary_terms,
+                                  temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
   return static_cast<int>(arguments.size());
 }
 
@@ -7206,51 +7161,32 @@ ExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 }
 
 void
-ExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                              bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                              const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                              const deriv_node_temp_terms_t &tef_terms) const
+ExternalFunctionNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                          const temporary_terms_t &temporary_terms,
+                                          const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                          const deriv_node_temp_terms_t &tef_terms) const
 {
-  if (auto this2 = const_cast<ExternalFunctionNode *>(this);
-      temporary_terms.contains(this2))
-    {
-      if (dynamic)
-        {
-          FLDT_ fldt(temporary_terms_idxs.at(this2));
-          fldt.write(CompileCode, instruction_number);
-        }
-      else
-        {
-          FLDST_ fldst(temporary_terms_idxs.at(this2));
-          fldst.write(CompileCode, instruction_number);
-        }
-      return;
-    }
+  if (checkIfTemporaryTermThenWriteBytecode(code_file, temporary_terms, temporary_terms_idxs, dynamic))
+    return;
 
   if (!lhs_rhs)
-    {
-      FLDTEF_ fldtef(getIndxInTefTerms(symb_id, tef_terms));
-      fldtef.write(CompileCode, instruction_number);
-    }
+    code_file << FLDTEF_{getIndxInTefTerms(symb_id, tef_terms)};
   else
-    {
-      FSTPTEF_ fstptef(getIndxInTefTerms(symb_id, tef_terms));
-      fstptef.write(CompileCode, instruction_number);
-    }
+    code_file << FSTPTEF_{getIndxInTefTerms(symb_id, tef_terms)};
 }
 
 void
-ExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                    const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                                    deriv_node_temp_terms_t &tef_terms) const
+ExternalFunctionNode::writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                                          const temporary_terms_t &temporary_terms,
+                                                          const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                          deriv_node_temp_terms_t &tef_terms) const
 {
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
   assert(first_deriv_symb_id != ExternalFunctionsTable::IDSetButNoNameProvided);
 
   for (auto argument : arguments)
-    argument->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                            temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+    argument->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs, temporary_terms,
+                                                  temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 
   if (!alreadyWrittenAsTefTerm(symb_id, tef_terms))
     {
@@ -7267,10 +7203,10 @@ ExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsign
         nb_output_arguments = 2;
       else
         nb_output_arguments = 1;
-      int nb_input_arguments{compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                                              temporary_terms_idxs, dynamic, steady_dynamic, tef_terms)};
+      int nb_input_arguments{writeBytecodeExternalFunctionArguments(code_file, lhs_rhs, temporary_terms,
+                                                                    temporary_terms_idxs, dynamic, steady_dynamic, tef_terms)};
 
-      FCALL_ fcall(nb_output_arguments, nb_input_arguments, datatree.symbol_table.getName(symb_id), indx);
+      FCALL_ fcall{nb_output_arguments, nb_input_arguments, datatree.symbol_table.getName(symb_id), indx};
       switch (nb_output_arguments)
         {
         case 1:
@@ -7283,9 +7219,7 @@ ExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsign
           fcall.set_function_type(ExternalFunctionType::withFirstAndSecondDerivative);
           break;
         }
-      fcall.write(CompileCode, instruction_number);
-      FSTPTEF_ fstptef(indx);
-      fstptef.write(CompileCode, instruction_number);
+      code_file << fcall << FSTPTEF_{indx};
     }
 }
 
@@ -7599,39 +7533,21 @@ FirstDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
 }
 
 void
-FirstDerivExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                        const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                        const deriv_node_temp_terms_t &tef_terms) const
+FirstDerivExternalFunctionNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                                    const temporary_terms_t &temporary_terms,
+                                                    const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                    const deriv_node_temp_terms_t &tef_terms) const
 {
-  if (auto this2 = const_cast<FirstDerivExternalFunctionNode *>(this);
-      temporary_terms.contains(this2))
-    {
-      if (dynamic)
-        {
-          FLDT_ fldt(temporary_terms_idxs.at(this2));
-          fldt.write(CompileCode, instruction_number);
-        }
-      else
-        {
-          FLDST_ fldst(temporary_terms_idxs.at(this2));
-          fldst.write(CompileCode, instruction_number);
-        }
-      return;
-    }
+  if (checkIfTemporaryTermThenWriteBytecode(code_file, temporary_terms, temporary_terms_idxs, dynamic))
+    return;
+
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
   assert(first_deriv_symb_id != ExternalFunctionsTable::IDSetButNoNameProvided);
 
   if (!lhs_rhs)
-    {
-      FLDTEFD_ fldtefd(getIndxInTefTerms(symb_id, tef_terms), inputIndex);
-      fldtefd.write(CompileCode, instruction_number);
-    }
+    code_file << FLDTEFD_{getIndxInTefTerms(symb_id, tef_terms), inputIndex};
   else
-    {
-      FSTPTEFD_ fstptefd(getIndxInTefTerms(symb_id, tef_terms), inputIndex);
-      fstptefd.write(CompileCode, instruction_number);
-    }
+    code_file << FSTPTEFD_{getIndxInTefTerms(symb_id, tef_terms), inputIndex};
 }
 
 void
@@ -7760,10 +7676,10 @@ FirstDerivExternalFunctionNode::writeJsonExternalFunctionOutput(vector<string> &
 }
 
 void
-FirstDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                                              bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                              const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                                              deriv_node_temp_terms_t &tef_terms) const
+FirstDerivExternalFunctionNode::writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                                                    const temporary_terms_t &temporary_terms,
+                                                                    const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                                    deriv_node_temp_terms_t &tef_terms) const
 {
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
   assert(first_deriv_symb_id != ExternalFunctionsTable::IDSetButNoNameProvided);
@@ -7773,30 +7689,28 @@ FirstDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCo
   if (first_deriv_symb_id == symb_id)
     {
       expr_t parent = datatree.AddExternalFunction(symb_id, arguments);
-      parent->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs,
-                                            temporary_terms, temporary_terms_idxs,
-                                            dynamic, steady_dynamic, tef_terms);
+      parent->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs,
+                                                  temporary_terms, temporary_terms_idxs,
+                                                  dynamic, steady_dynamic, tef_terms);
       return;
     }
 
   if (alreadyWrittenAsTefTerm(first_deriv_symb_id, tef_terms))
     return;
 
-  int nb_add_input_arguments{compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                                              temporary_terms_idxs, dynamic, steady_dynamic, tef_terms)};
+  int nb_add_input_arguments{writeBytecodeExternalFunctionArguments(code_file, lhs_rhs, temporary_terms,
+                                                                    temporary_terms_idxs, dynamic, steady_dynamic, tef_terms)};
   if (first_deriv_symb_id == ExternalFunctionsTable::IDNotSet)
     {
       int nb_input_arguments{0};
       int nb_output_arguments{1};
       int indx = getIndxInTefTerms(symb_id, tef_terms);
-      FCALL_ fcall(nb_output_arguments, nb_input_arguments, "jacob_element", indx);
+      FCALL_ fcall{nb_output_arguments, nb_input_arguments, "jacob_element", indx};
       fcall.set_arg_func_name(datatree.symbol_table.getName(symb_id));
       fcall.set_row(inputIndex);
       fcall.set_nb_add_input_arguments(nb_add_input_arguments);
       fcall.set_function_type(ExternalFunctionType::numericalFirstDerivative);
-      fcall.write(CompileCode, instruction_number);
-      FSTPTEFD_ fstptefd(indx, inputIndex);
-      fstptefd.write(CompileCode, instruction_number);
+      code_file << fcall << FSTPTEFD_{indx, inputIndex};
     }
   else
     {
@@ -7807,11 +7721,9 @@ FirstDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCo
 
       int nb_output_arguments{1};
 
-      FCALL_ fcall(nb_output_arguments, nb_add_input_arguments, datatree.symbol_table.getName(first_deriv_symb_id), indx);
+      FCALL_ fcall{nb_output_arguments, nb_add_input_arguments, datatree.symbol_table.getName(first_deriv_symb_id), indx};
       fcall.set_function_type(ExternalFunctionType::firstDerivative);
-      fcall.write(CompileCode, instruction_number);
-      FSTPTEFD_ fstptefd(indx, inputIndex);
-      fstptefd.write(CompileCode, instruction_number);
+      code_file << fcall << FSTPTEFD_{indx, inputIndex};
     }
 }
 
@@ -8138,47 +8050,28 @@ SecondDerivExternalFunctionNode::computeXrefs(EquationInfo &ei) const
 }
 
 void
-SecondDerivExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                                         bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                         const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                         const deriv_node_temp_terms_t &tef_terms) const
+SecondDerivExternalFunctionNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                                     const temporary_terms_t &temporary_terms,
+                                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                     const deriv_node_temp_terms_t &tef_terms) const
 {
-  if (auto this2 = const_cast<SecondDerivExternalFunctionNode *>(this);
-      temporary_terms.contains(this2))
-    {
-      if (dynamic)
-        {
-          FLDT_ fldt(temporary_terms_idxs.at(this2));
-          fldt.write(CompileCode, instruction_number);
-        }
-      else
-        {
-          FLDST_ fldst(temporary_terms_idxs.at(this2));
-          fldst.write(CompileCode, instruction_number);
-        }
-      return;
-    }
+  if (checkIfTemporaryTermThenWriteBytecode(code_file, temporary_terms, temporary_terms_idxs, dynamic))
+    return;
 
   int second_deriv_symb_id = datatree.external_functions_table.getSecondDerivSymbID(symb_id);
   assert(second_deriv_symb_id != ExternalFunctionsTable::IDSetButNoNameProvided);
 
   if (!lhs_rhs)
-    {
-      FLDTEFDD_ fldtefdd(getIndxInTefTerms(symb_id, tef_terms), inputIndex1, inputIndex2);
-      fldtefdd.write(CompileCode, instruction_number);
-    }
+    code_file << FLDTEFDD_{getIndxInTefTerms(symb_id, tef_terms), inputIndex1, inputIndex2};
   else
-    {
-      FSTPTEFDD_ fstptefdd(getIndxInTefTerms(symb_id, tef_terms), inputIndex1, inputIndex2);
-      fstptefdd.write(CompileCode, instruction_number);
-    }
+    code_file << FSTPTEFDD_{getIndxInTefTerms(symb_id, tef_terms), inputIndex1, inputIndex2};
 }
 
 void
-SecondDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                                               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                                               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                                               deriv_node_temp_terms_t &tef_terms) const
+SecondDerivExternalFunctionNode::writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                                                     const temporary_terms_t &temporary_terms,
+                                                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                                     deriv_node_temp_terms_t &tef_terms) const
 {
   int second_deriv_symb_id = datatree.external_functions_table.getSecondDerivSymbID(symb_id);
   assert(second_deriv_symb_id != ExternalFunctionsTable::IDSetButNoNameProvided);
@@ -8188,31 +8081,29 @@ SecondDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileC
   if (second_deriv_symb_id == symb_id)
     {
       expr_t parent = datatree.AddExternalFunction(symb_id, arguments);
-      parent->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs,
-                                            temporary_terms, temporary_terms_idxs,
-                                            dynamic, steady_dynamic, tef_terms);
+      parent->writeBytecodeExternalFunctionOutput(code_file, lhs_rhs,
+                                                  temporary_terms, temporary_terms_idxs,
+                                                  dynamic, steady_dynamic, tef_terms);
       return;
     }
 
   if (alreadyWrittenAsTefTerm(second_deriv_symb_id, tef_terms))
     return;
 
-  int nb_add_input_arguments{compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs, temporary_terms,
-                                                              temporary_terms_idxs, dynamic, steady_dynamic, tef_terms)};
+  int nb_add_input_arguments{writeBytecodeExternalFunctionArguments(code_file, lhs_rhs, temporary_terms,
+                                                                    temporary_terms_idxs, dynamic, steady_dynamic, tef_terms)};
   if (second_deriv_symb_id == ExternalFunctionsTable::IDNotSet)
     {
       int nb_input_arguments{0};
       int nb_output_arguments{1};
       int indx = getIndxInTefTerms(symb_id, tef_terms);
-      FCALL_ fcall(nb_output_arguments, nb_input_arguments, "hess_element", indx);
+      FCALL_ fcall{nb_output_arguments, nb_input_arguments, "hess_element", indx};
       fcall.set_arg_func_name(datatree.symbol_table.getName(symb_id));
       fcall.set_row(inputIndex1);
       fcall.set_col(inputIndex2);
       fcall.set_nb_add_input_arguments(nb_add_input_arguments);
       fcall.set_function_type(ExternalFunctionType::numericalSecondDerivative);
-      fcall.write(CompileCode, instruction_number);
-      FSTPTEFDD_ fstptefdd(indx, inputIndex1, inputIndex2);
-      fstptefdd.write(CompileCode, instruction_number);
+      code_file << fcall << FSTPTEFDD_{indx, inputIndex1, inputIndex2};
     }
   else
     {
@@ -8221,11 +8112,9 @@ SecondDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileC
 
       int nb_output_arguments{1};
 
-      FCALL_ fcall(nb_output_arguments, nb_add_input_arguments, datatree.symbol_table.getName(second_deriv_symb_id), indx);
+      FCALL_ fcall{nb_output_arguments, nb_add_input_arguments, datatree.symbol_table.getName(second_deriv_symb_id), indx};
       fcall.set_function_type(ExternalFunctionType::secondDerivative);
-      fcall.write(CompileCode, instruction_number);
-      FSTPTEFDD_ fstptefdd(indx, inputIndex1, inputIndex2);
-      fstptefdd.write(CompileCode, instruction_number);
+      code_file << fcall << FSTPTEFDD_{indx, inputIndex1, inputIndex2};
     }
 }
 
@@ -8428,10 +8317,10 @@ SubModelNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &
 }
 
 void
-SubModelNode::compile(ostream &CompileCode, unsigned int &instruction_number,
-                            bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                            const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                            const deriv_node_temp_terms_t &tef_terms) const
+SubModelNode::writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs,
+                                  const temporary_terms_t &temporary_terms,
+                                  const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                  const deriv_node_temp_terms_t &tef_terms) const
 {
   cerr << "SubModelNode::compile not implemented." << endl;
   exit(EXIT_FAILURE);
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 97d58d4cf17bcdfe81ef5d1bcb227b13cc9ad84d..b08320eab0b0168b6b849924b93a6ca07af11cc9 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -31,6 +31,7 @@ using namespace std;
 
 #include "CommonEnums.hh"
 #include "ExternalFunctionsTable.hh"
+#include "Bytecode.hh"
 
 class DataTree;
 class NumConstNode;
@@ -240,6 +241,12 @@ protected:
                                      const temporary_terms_t &temporary_terms,
                                      const temporary_terms_idxs_t &temporary_terms_idxs) const;
 
+  // Same as above, for the bytecode case
+  bool checkIfTemporaryTermThenWriteBytecode(BytecodeWriter &code_file,
+                                             const temporary_terms_t &temporary_terms,
+                                             const temporary_terms_idxs_t &temporary_terms_idxs,
+                                             bool dynamic) const;
+
   // Internal helper for matchVariableTimesConstantTimesParam()
   virtual void matchVTCTPHelper(optional<int> &var_id, int &lag, optional<int> &param_id, double &constant, bool at_denominator) const;
 
@@ -360,10 +367,10 @@ public:
                                                deriv_node_temp_terms_t &tef_terms,
                                                bool isdynamic = true) const;
 
-  virtual void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                             bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                             const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                             deriv_node_temp_terms_t &tef_terms) const;
+  virtual void writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file,
+                                                   bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                   const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                                   deriv_node_temp_terms_t &tef_terms) const;
 
   //! Computes the set of all variables of a given symbol type in the expression (with information on lags)
   /*!
@@ -406,7 +413,7 @@ public:
   };
 
   virtual double eval(const eval_context_t &eval_context) const noexcept(false) = 0;
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const = 0;
+  virtual void writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const = 0;
 
   //! Creates a static version of this node
   /*!
@@ -788,7 +795,7 @@ public:
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
@@ -860,7 +867,7 @@ public:
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   SymbolType get_type() const;
@@ -955,7 +962,7 @@ public:
                                        const temporary_terms_t &temporary_terms,
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
-  void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+  void writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file,
                                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
                                      const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
                                      deriv_node_temp_terms_t &tef_terms) const override;
@@ -963,7 +970,7 @@ public:
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   static double eval_opcode(UnaryOpcode op_code, double v) noexcept(false);
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
@@ -1058,15 +1065,15 @@ public:
                                        const temporary_terms_t &temporary_terms,
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
-  void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                     deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file,
+                                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                           deriv_node_temp_terms_t &tef_terms) const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   static double eval_opcode(double v1, BinaryOpcode op_code, double v2, int derivOrder) noexcept(false);
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
@@ -1199,15 +1206,15 @@ public:
                                        const temporary_terms_t &temporary_terms,
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
-  void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                     deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file,
+                                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                           deriv_node_temp_terms_t &tef_terms) const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) noexcept(false);
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
@@ -1281,6 +1288,10 @@ protected:
   void writeExternalFunctionArguments(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
   void writeJsonASTExternalFunctionArguments(ostream &output) const;
   void writeJsonExternalFunctionArguments(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const;
+  int writeBytecodeExternalFunctionArguments(BytecodeWriter &code_file,
+                                             bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                             const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                             const deriv_node_temp_terms_t &tef_terms) const;
   /*! Returns a predicate that tests whether an other ExprNode is an external
     function which is computed by the same external function call (i.e. it has
     the same so-called "Tef" index) */
@@ -1307,19 +1318,14 @@ public:
                                        const temporary_terms_t &temporary_terms,
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic = true) const override = 0;
-  void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                     deriv_node_temp_terms_t &tef_terms) const override = 0;
+  void writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file,
+                                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                           deriv_node_temp_terms_t &tef_terms) const override = 0;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   double eval(const eval_context_t &eval_context) const noexcept(false) override;
-  int compileExternalFunctionArguments(ostream &CompileCode, unsigned int &instruction_number,
-                                       bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                       const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                       const deriv_node_temp_terms_t &tef_terms) const;
-
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override = 0;
+  void writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override = 0;
   expr_t toStatic(DataTree &static_datatree) const override = 0;
   void computeXrefs(EquationInfo &ei) const override = 0;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
@@ -1393,11 +1399,11 @@ public:
                                        const temporary_terms_t &temporary_terms,
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
-  void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                     deriv_node_temp_terms_t &tef_terms) const override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file,
+                                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                           deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file, bool lhs_rhs, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic, const deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const override;
@@ -1420,10 +1426,10 @@ public:
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number,
-               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-               const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file,
+                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                           const deriv_node_temp_terms_t &tef_terms) const override;
   void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                    const temporary_terms_t &temporary_terms,
                                    const temporary_terms_idxs_t &temporary_terms_idxs,
@@ -1432,10 +1438,10 @@ public:
                                        const temporary_terms_t &temporary_terms,
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
-  void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                     deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file,
+                                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                           deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const override;
@@ -1460,10 +1466,10 @@ public:
   void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
   void writeJsonAST(ostream &output) const override;
   void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number,
-               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-               const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file,
+                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                           const deriv_node_temp_terms_t &tef_terms) const override;
   void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                    const temporary_terms_t &temporary_terms,
                                    const temporary_terms_idxs_t &temporary_terms_idxs,
@@ -1472,10 +1478,10 @@ public:
                                        const temporary_terms_t &temporary_terms,
                                        deriv_node_temp_terms_t &tef_terms,
                                        bool isdynamic) const override;
-  void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
-                                     bool lhs_rhs, const temporary_terms_t &temporary_terms,
-                                     const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-                                     deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeExternalFunctionOutput(BytecodeWriter &code_file,
+                                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                                           deriv_node_temp_terms_t &tef_terms) const override;
   expr_t toStatic(DataTree &static_datatree) const override;
   void computeXrefs(EquationInfo &ei) const override;
   expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const override;
@@ -1528,10 +1534,10 @@ public:
   expr_t substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   void computeSubExprContainingVariable(int symb_id, int lag, set<expr_t> &contain_var) const override;
   BinaryOpNode *normalizeEquationHelper(const set<expr_t> &contain_var, expr_t rhs) const override;
-  void compile(ostream &CompileCode, unsigned int &instruction_number,
-               bool lhs_rhs, const temporary_terms_t &temporary_terms,
-               const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
-               const deriv_node_temp_terms_t &tef_terms) const override;
+  void writeBytecodeOutput(BytecodeWriter &code_file,
+                           bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                           const temporary_terms_idxs_t &temporary_terms_idxs, bool dynamic, bool steady_dynamic,
+                           const deriv_node_temp_terms_t &tef_terms) const override;
   void collectVARLHSVariable(set<expr_t> &result) const override;
   void collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &result) const override;
   bool isNumConstNodeEqualTo(double value) const override;
diff --git a/src/Makefile.am b/src/Makefile.am
index d3b4bcc2665095702c9637148700bae0c9043929..45d6d7e54dd26e1e990d300e41112410761541eb 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -46,6 +46,7 @@ dynare_preprocessor_SOURCES = \
 	VariableDependencyGraph.hh \
 	DynareMain.cc \
 	MacroExpandModFile.cc \
+	Bytecode.cc \
 	Bytecode.hh \
 	ExternalFunctionsTable.cc \
 	ExternalFunctionsTable.hh \
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index ee58009be311687a6fafed3d45695d26df4ed5ab..0e0e74f9b455a0519bb60cdc7a8e563648f053be 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -19,7 +19,6 @@
 
 #include "ModelTree.hh"
 #include "VariableDependencyGraph.hh"
-#include "Bytecode.hh"
 
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wold-style-cast"
@@ -1256,27 +1255,20 @@ ModelTree::testNestedParenthesis(const string &str) const
 }
 
 void
-ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, bool dynamic, bool steady_dynamic, temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const
+ModelTree::writeBytecodeTemporaryTerms(BytecodeWriter &code_file, bool dynamic, bool steady_dynamic, 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->compileExternalFunctionOutput(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+        tt->writeBytecodeExternalFunctionOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 
-      FNUMEXPR_ fnumexpr(ExpressionType::TemporaryTerm, idx);
-      fnumexpr.write(code_file, instruction_number);
-      tt->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+      code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, idx};
+      tt->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
       if (dynamic)
-        {
-          FSTPT_ fstpt(idx);
-          fstpt.write(code_file, instruction_number);
-        }
+        code_file << FSTPT_{idx};
       else
-        {
-          FSTPST_ fstpst(idx);
-          fstpst.write(code_file, instruction_number);
-        }
+        code_file << FSTPST_{idx};
     }
 }
 
@@ -1392,15 +1384,14 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type,
 }
 
 void
-ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_number, bool dynamic, bool steady_dynamic, const temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
+ModelTree::writeBytecodeModelEquations(BytecodeWriter &code_file, bool dynamic, bool steady_dynamic, 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;
       expr_t rhs = eq_node->arg2;
-      FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, eq);
-      fnumexpr.write(code_file, instruction_number);
+      code_file << FNUMEXPR_{ExpressionType::ModelEquation, eq};
       // Test if the right hand side of the equation is empty.
       double vrhs = 1.0;
       try
@@ -1413,20 +1404,15 @@ ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_n
 
       if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
         {
-          lhs->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-          rhs->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-
-          FBINARY_ fbinary{BinaryOpcode::minus};
-          fbinary.write(code_file, instruction_number);
+          lhs->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+          rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
 
-          FSTPR_ fstpr(eq);
-          fstpr.write(code_file, instruction_number);
+          code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{eq};
         }
       else // The right hand side of the equation is empty ==> residual=lhs;
         {
-          lhs->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
-          FSTPR_ fstpr(eq);
-          fstpr.write(code_file, instruction_number);
+          lhs->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, dynamic, steady_dynamic, tef_terms);
+          code_file << FSTPR_{eq};
         }
     }
 }
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 32c119745bc84618a22127bffe4c42e3f5a94e43..63c3491e1ec9f86b3554965cbad81bf82323d4d6 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -32,6 +32,7 @@
 #include "DataTree.hh"
 #include "EquationTags.hh"
 #include "ExtendedPreprocessorTypes.hh"
+#include "Bytecode.hh"
 
 using namespace std;
 
@@ -236,8 +237,8 @@ protected:
   //! Writes temporary terms
   void writeTemporaryTerms(const temporary_terms_t &tt, temporary_terms_t &temp_term_union, const temporary_terms_idxs_t &tt_idxs, ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
   void writeJsonTemporaryTerms(const temporary_terms_t &tt, temporary_terms_t &temp_term_union, ostream &output, deriv_node_temp_terms_t &tef_terms, const string &concat) const;
-  //! Compiles temporary terms
-  void compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, bool dynamic, bool steady_dynamic, temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
+  //! Writes temporary terms in bytecode
+  void writeBytecodeTemporaryTerms(BytecodeWriter &code_file, bool dynamic, bool steady_dynamic, 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;
   //! Fixes output when there are more than 32 nested parens, Issue #1201
@@ -259,8 +260,8 @@ 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;
-  //! Compiles model equations
-  void compileModelEquations(ostream &code_file, unsigned int &instruction_number, bool dynamic, bool steady_dynamic, const temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
+  //! Writes model equations in bytecode
+  void writeBytecodeModelEquations(BytecodeWriter &code_file, bool dynamic, bool steady_dynamic, const temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, 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;
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 6ae1b984b02c1a0a7a5b5aee8a9b9c26ca2f576e..13a4107ef77a996a70276fcbc80fc30aaab48974 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -26,7 +26,6 @@
 
 #include "StaticModel.hh"
 #include "DynamicModel.hh"
-#include "Bytecode.hh"
 
 void
 StaticModel::copyHelper(const StaticModel &m)
@@ -104,29 +103,23 @@ StaticModel::StaticModel(const DynamicModel &m) :
 }
 
 void
-StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, 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
+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->compile(code_file, instruction_number, false, temporary_terms, temporary_terms_idxs, false, false, tef_terms);
+    it->second->writeBytecodeOutput(code_file, false, temporary_terms, temporary_terms_idxs, false, false, tef_terms);
   else
-    {
-      FLDZ_ fldz;
-      fldz.write(code_file, instruction_number);
-    }
+    code_file << FLDZ_{};
 }
 
 void
-StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, 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
+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->compile(code_file, instruction_number, false, temporary_terms, temporary_terms_idxs, false, false, tef_terms);
+    it->second->writeBytecodeOutput(code_file, false, temporary_terms, temporary_terms_idxs, false, false, tef_terms);
   else
-    {
-      FLDZ_ fldz;
-      fldz.write(code_file, instruction_number);
-    }
+    code_file << FLDZ_{};
 }
 
 void
@@ -381,17 +374,9 @@ void
 StaticModel::writeStaticBytecode(const string &basename) const
 {
   ostringstream tmp_output;
-  ofstream code_file;
-  unsigned int instruction_number = 0;
   bool file_open = false;
 
-  string main_name = basename + "/model/bytecode/static.cod";
-  code_file.open(main_name, ios::out | ios::binary | ios::ate);
-  if (!code_file.is_open())
-    {
-      cerr << R"(Error : Can't open file ")" << main_name << R"(" for writing)" << endl;
-      exit(EXIT_FAILURE);
-    }
+  BytecodeWriter code_file{basename + "/model/bytecode/static.cod"};
   int count_u;
   int u_count_int = 0;
 
@@ -403,37 +388,32 @@ StaticModel::writeStaticBytecode(const string &basename) const
   temporary_terms.insert(temporary_terms_derivatives[1].begin(), temporary_terms_derivatives[1].end());
 
   //Temporary variables declaration
-  FDIMST_ fdimst{static_cast<int>(temporary_terms.size())};
-  fdimst.write(code_file, instruction_number);
-  FBEGINBLOCK_ 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());
-  fbeginblock.write(code_file, instruction_number);
+  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;
 
-  compileTemporaryTerms(code_file, instruction_number, false, false, temporary_terms_union, temporary_terms_idxs, tef_terms);
+  writeBytecodeTemporaryTerms(code_file, false, false, temporary_terms_union, temporary_terms_idxs, tef_terms);
 
-  compileModelEquations(code_file, instruction_number, false, false, temporary_terms_union, temporary_terms_idxs, tef_terms);
+  writeBytecodeModelEquations(code_file, false, false, temporary_terms_union, temporary_terms_idxs, tef_terms);
 
-  FENDEQU_ fendequ;
-  fendequ.write(code_file, instruction_number);
+  code_file << FENDEQU_{};
 
-  // Get the current code_file position and jump if eval = true
-  streampos pos1 = code_file.tellp();
-  FJMPIFEVAL_ fjmp_if_eval(0);
-  fjmp_if_eval.write(code_file, instruction_number);
-  int prev_instruction_number = instruction_number;
+  // 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();
@@ -445,57 +425,41 @@ StaticModel::writeStaticBytecode(const string &basename) const
           int eq = indices[0];
           int symb = getSymbIDByDerivID(deriv_id);
           int var = symbol_table.getTypeSpecificID(symb);
-          FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eq, var);
-          fnumexpr.write(code_file, instruction_number);
+          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->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, false, false, tef_terms);
+          d1->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, false, false, tef_terms);
 
-          FSTPSU_ fstpsu(count_u);
-          fstpsu.write(code_file, instruction_number);
+          code_file << FSTPSU_{count_u};
           count_u++;
         }
     }
   for (int i = 0; i < symbol_table.endo_nbr(); i++)
     {
-      FLDR_ fldr(i);
-      fldr.write(code_file, instruction_number);
+      code_file << FLDR_{i};
       if (my_derivatives[i].size())
         {
           for (bool printed_something{false};
                const auto &it : my_derivatives[i])
             {
-              FLDSU_ fldsu(it.second);
-              fldsu.write(code_file, instruction_number);
-              FLDSV_ fldsv{SymbolType::endogenous, it.first};
-              fldsv.write(code_file, instruction_number);
-              FBINARY_ fbinary{BinaryOpcode::times};
-              fbinary.write(code_file, instruction_number);
+              code_file << FLDSU_{it.second}
+                << FLDSV_{SymbolType::endogenous, it.first}
+                << FBINARY_{BinaryOpcode::times};
               if (exchange(printed_something, true))
-                {
-                  FBINARY_ fbinary{BinaryOpcode::plus};
-                  fbinary.write(code_file, instruction_number);
-                }
+                code_file << FBINARY_{BinaryOpcode::plus};
             }
-          FBINARY_ fbinary{BinaryOpcode::minus};
-          fbinary.write(code_file, instruction_number);
+          code_file << FBINARY_{BinaryOpcode::minus};
         }
-      FSTPSU_ fstpsu(i);
-      fstpsu.write(code_file, instruction_number);
+      code_file << FSTPSU_{i};
     }
-  // Get the current code_file position and jump = true
-  streampos pos2 = code_file.tellp();
-  FJMP_ fjmp(0);
-  fjmp.write(code_file, instruction_number);
-  // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
-  streampos pos3 = code_file.tellp();
-  code_file.seekp(pos1);
-  FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
-  fjmp_if_eval1.write(code_file, instruction_number);
-  code_file.seekp(pos3);
-  prev_instruction_number = instruction_number;
+
+  // 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;
 
@@ -508,30 +472,21 @@ StaticModel::writeStaticBytecode(const string &basename) const
           int eq = indices[0];
           int symb = getSymbIDByDerivID(deriv_id);
           int var = symbol_table.getTypeSpecificID(symb);
-          FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eq, var);
-          fnumexpr.write(code_file, instruction_number);
+          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->compile(code_file, instruction_number, false, temporary_terms_union, temporary_terms_idxs, false, false, tef_terms);
-          FSTPG2_ fstpg2(eq, var);
-          fstpg2.write(code_file, instruction_number);
+          d1->writeBytecodeOutput(code_file, false, temporary_terms_union, temporary_terms_idxs, false, false, tef_terms);
+          code_file << FSTPG2_{eq, var};
         }
     }
 
-  // Set codefile position to previous JMP_ and set the number of instructions to jump
-  pos1 = code_file.tellp();
-  code_file.seekp(pos2);
-  FJMP_ fjmp1(instruction_number - prev_instruction_number);
-  fjmp1.write(code_file, instruction_number);
-  code_file.seekp(pos1);
-
-  FENDBLOCK_ fendblock;
-  fendblock.write(code_file, instruction_number);
-  FEND_ fend;
-  fend.write(code_file, instruction_number);
-  code_file.close();
+  // 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_{};
 }
 
 void
@@ -551,8 +506,6 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
   int i, v;
   string tmp_s;
   ostringstream tmp_output;
-  ofstream code_file;
-  unsigned int instruction_number = 0;
   expr_t lhs = nullptr, rhs = nullptr;
   BinaryOpNode *eq_node;
   Uff Uf[symbol_table.endo_nbr()];
@@ -560,17 +513,10 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
   vector<int> feedback_variables;
   bool file_open = false;
 
-  string main_name = basename + "/model/bytecode/static.cod";
-  code_file.open(main_name, ios::out | ios::binary | ios::ate);
-  if (!code_file.is_open())
-    {
-      cerr << R"(Error : Can't open file ")" << main_name << R"(" for writing)" << endl;
-      exit(EXIT_FAILURE);
-    }
-  //Temporary variables declaration
+  BytecodeWriter code_file{basename + "/model/bytecode/static.cod"};
 
-  FDIMST_ fdimst{static_cast<int>(blocks_temporary_terms_idxs.size())};
-  fdimst.write(code_file, instruction_number);
+  //Temporary variables declaration
+  code_file << FDIMST_{static_cast<int>(blocks_temporary_terms_idxs.size())};
 
   temporary_terms_t temporary_terms_union;
 
@@ -578,10 +524,7 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
     {
       feedback_variables.clear();
       if (block > 0)
-        {
-          FENDBLOCK_ fendblock;
-          fendblock.write(code_file, instruction_number);
-        }
+        code_file << FENDBLOCK_{};
       int count_u;
       int u_count_int = 0;
       BlockSimulationType simulation_type = blocks[block].simulation_type;
@@ -598,26 +541,22 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
           file_open = true;
         }
 
-      FBEGINBLOCK_ 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);
-
-      fbeginblock.write(code_file, instruction_number);
-
-      // Get the current code_file position and jump if eval = true
-      streampos pos1 = code_file.tellp();
-      FJMPIFEVAL_ fjmp_if_eval(0);
-      fjmp_if_eval.write(code_file, instruction_number);
-      int prev_instruction_number = instruction_number;
+      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;
@@ -630,13 +569,11 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
                            for (auto it : blocks_temporary_terms[block][eq])
                              {
                                if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-                                 it->compileExternalFunctionOutput(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                                 it->writeBytecodeExternalFunctionOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
 
-                               FNUMEXPR_ fnumexpr{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)};
-                               fnumexpr.write(code_file, instruction_number);
-                               it->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-                               FSTPST_ fstpst{blocks_temporary_terms_idxs.at(it)};
-                               fstpst.write(code_file, instruction_number);
+                               code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)};
+                               it->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                               code_file << FSTPST_{blocks_temporary_terms_idxs.at(it)};
                                temporary_terms_union.insert(it);
                              }
                          };
@@ -654,25 +591,22 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
             case BlockSimulationType::evaluateBackward:
             case BlockSimulationType::evaluateForward:
               equ_type = getBlockEquationType(block, i);
-              {
-                FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, getBlockEquationID(block, i));
-                fnumexpr.write(code_file, instruction_number);
-              }
+              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->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                  rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                  lhs->writeBytecodeOutput(code_file, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
                 }
               else if (equ_type == EquationType::evaluateRenormalized)
                 {
                   eq_node = getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                  rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                  lhs->writeBytecodeOutput(code_file, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
                 }
               break;
             case BlockSimulationType::solveBackwardComplete:
@@ -686,23 +620,17 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
               goto end;
             default:
             end:
-              FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, getBlockEquationID(block, i));
-              fnumexpr.write(code_file, instruction_number);
+              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
               eq_node = getBlockEquationExpr(block, i);
               lhs = eq_node->arg1;
               rhs = eq_node->arg2;
-              lhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-              rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-
-              FBINARY_ fbinary{BinaryOpcode::minus};
-              fbinary.write(code_file, instruction_number);
+              lhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+              rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
 
-              FSTPR_ fstpr(i - block_recursive);
-              fstpr.write(code_file, instruction_number);
+              code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
             }
         }
-      FENDEQU_ fendequ;
-      fendequ.write(code_file, instruction_number);
+      code_file << FENDEQU_{};
 
       // The Jacobian if we have to solve the block
       if (simulation_type != BlockSimulationType::evaluateBackward
@@ -715,15 +643,9 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
             {
             case BlockSimulationType::solveBackwardSimple:
             case BlockSimulationType::solveForwardSimple:
-              {
-                FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, 0, 0);
-                fnumexpr.write(code_file, instruction_number);
-              }
-              compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-              {
-                FSTPG_ fstpg(0);
-                fstpg.write(code_file, instruction_number);
-              }
+              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:
@@ -749,37 +671,23 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
                       Uf[eqr].Ufl->pNext = nullptr;
                       Uf[eqr].Ufl->u = count_u;
                       Uf[eqr].Ufl->var = varr;
-                      FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eqr, varr);
-                      fnumexpr.write(code_file, instruction_number);
-                      compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-                      FSTPSU_ fstpsu(count_u);
-                      fstpsu.write(code_file, instruction_number);
+                      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)
                   {
-                    FLDR_ fldr(i-block_recursive);
-                    fldr.write(code_file, instruction_number);
-
-                    FLDZ_ fldz;
-                    fldz.write(code_file, instruction_number);
+                    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)
-                      {
-                        FLDSU_ fldsu(Uf[v].Ufl->u);
-                        fldsu.write(code_file, instruction_number);
-                        FLDSV_ fldsv{SymbolType::endogenous, Uf[v].Ufl->var};
-                        fldsv.write(code_file, instruction_number);
-
-                        FBINARY_ fbinary{BinaryOpcode::times};
-                        fbinary.write(code_file, instruction_number);
-
-                        FCUML_ fcuml;
-                        fcuml.write(code_file, instruction_number);
-                      }
+                      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)
                       {
@@ -787,12 +695,8 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
                         free(Uf[v].Ufl);
                         Uf[v].Ufl = Uf[v].Ufl_First;
                       }
-                    FBINARY_ fbinary{BinaryOpcode::minus};
-                    fbinary.write(code_file, instruction_number);
-
-                    FSTPSU_ fstpsu(i - block_recursive);
-                    fstpsu.write(code_file, instruction_number);
-
+                    code_file << FBINARY_{BinaryOpcode::minus}
+                      << FSTPSU_{i - block_recursive};
                   }
               break;
             default:
@@ -800,17 +704,11 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
             }
         }
 
-      // Get the current code_file position and jump = true
-      streampos pos2 = code_file.tellp();
-      FJMP_ fjmp(0);
-      fjmp.write(code_file, instruction_number);
-      // Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
-      streampos pos3 = code_file.tellp();
-      code_file.seekp(pos1);
-      FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
-      fjmp_if_eval1.write(code_file, instruction_number);
-      code_file.seekp(pos3);
-      prev_instruction_number = instruction_number;
+      // 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;
@@ -828,25 +726,22 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
             case BlockSimulationType::evaluateBackward:
             case BlockSimulationType::evaluateForward:
               equ_type = getBlockEquationType(block, i);
-              {
-                FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, getBlockEquationID(block, i));
-                fnumexpr.write(code_file, instruction_number);
-              }
+              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->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                  rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                  lhs->writeBytecodeOutput(code_file, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
                 }
               else if (equ_type == EquationType::evaluateRenormalized)
                 {
                   eq_node = getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->arg1;
                   rhs = eq_node->arg2;
-                  rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-                  lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                  rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+                  lhs->writeBytecodeOutput(code_file, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
                 }
               break;
             case BlockSimulationType::solveBackwardComplete:
@@ -860,23 +755,17 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
               goto end_l;
             default:
             end_l:
-              FNUMEXPR_ fnumexpr(ExpressionType::ModelEquation, getBlockEquationID(block, i));
-              fnumexpr.write(code_file, instruction_number);
+              code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)};
               eq_node = getBlockEquationExpr(block, i);
               lhs = eq_node->arg1;
               rhs = eq_node->arg2;
-              lhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-              rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+              lhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+              rhs->writeBytecodeOutput(code_file, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
 
-              FBINARY_ fbinary{BinaryOpcode::minus};
-              fbinary.write(code_file, instruction_number);
-
-              FSTPR_ fstpr(i - block_recursive);
-              fstpr.write(code_file, instruction_number);
+              code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive};
             }
         }
-      FENDEQU_ fendequ_l;
-      fendequ_l.write(code_file, instruction_number);
+      code_file << FENDEQU_{};
 
       // The Jacobian if we have to solve the block determinsitic bloc
 
@@ -887,15 +776,9 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
         {
         case BlockSimulationType::solveBackwardSimple:
         case BlockSimulationType::solveForwardSimple:
-          {
-            FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, 0, 0);
-            fnumexpr.write(code_file, instruction_number);
-          }
-          compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-          {
-            FSTPG2_ fstpg2(0, 0);
-            fstpg2.write(code_file, instruction_number);
-          }
+          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:
@@ -907,30 +790,20 @@ StaticModel::writeStaticBlockBytecode(const string &basename) const
               auto &[eq, var, ignore] = indices;
               int eqr = getBlockEquationID(block, eq);
               int varr = getBlockVariableID(block, var);
-              FNUMEXPR_ fnumexpr(ExpressionType::FirstEndoDerivative, eqr, varr, 0);
-              fnumexpr.write(code_file, instruction_number);
-
-              compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms);
-
-              FSTPG2_ fstpg2(eq, var);
-              fstpg2.write(code_file, instruction_number);
+              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
-      pos1 = code_file.tellp();
-      code_file.seekp(pos2);
-      FJMP_ fjmp1(instruction_number - prev_instruction_number);
-      fjmp1.write(code_file, instruction_number);
-      code_file.seekp(pos1);
+      // 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});
     }
-  FENDBLOCK_ fendblock;
-  fendblock.write(code_file, instruction_number);
-  FEND_ fend;
-  fend.write(code_file, instruction_number);
-  code_file.close();
+  code_file << FENDBLOCK_{} << FEND_{};
 }
 
 void
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index b582516a313fd00f0d3840f63cd3e95bb38e5a79..6966c6f2b20d473fbe7c514dea07d0672cb9332c 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -24,6 +24,7 @@
 #include <filesystem>
 
 #include "ModelTree.hh"
+#include "Bytecode.hh"
 
 using namespace std;
 
@@ -77,10 +78,10 @@ private:
   */
   void evaluateJacobian(const eval_context_t &eval_context, jacob_map_t *j_m, bool dynamic);
 
-  //! Write derivative code of an equation w.r. to a variable
-  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, 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 code of an equation w.r. to a variable
-  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const;
+  //! 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;
 
   //! Get the type corresponding to a derivation ID
   SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;