diff --git a/CodeInterpreter.hh b/CodeInterpreter.hh
index b8b65a5300050bb2391326ca1e0f2302189c2fcd..3f1b55523ff08e7e411f3799c04e0193457371d7 100644
--- a/CodeInterpreter.hh
+++ b/CodeInterpreter.hh
@@ -93,7 +93,17 @@ enum Tags
 
     FOK,          //!< Used for debugging purpose - 21 (33)
 
-    FNUMEXPR      //!< Store the expression type and references - 22 (34)
+    FNUMEXPR,     //!< Store the expression type and references - 22 (34)
+
+    FCALL,        //!< Call an external function - 23 (35)
+    FPUSH,        //!< Push a double in the stack - 24 (36)
+    FPOP,         //!< Pop a double from the stack - 25 (37)
+    FLDTEF,       //!< Stores the result of an external function in the stack - 26 (38)
+    FSTPTEF,      //!< Loads the result of an external function from the stack- 27 (39)
+    FLDTEFD,      //!< Stores the result of an external function in the stack - 28 (40)
+    FSTPTEFD,     //!< Loads the result of an external function from the stack- 29 (41)
+    FLDTEFDD,     //!< Stores the result of an external function in the stack - 28 (42)
+    FSTPTEFDD     //!< Loads the result of an external function from the stack- 29 (43)
 
   };
 
@@ -207,6 +217,17 @@ enum TrinaryOpcode
     oNormpdf
   };
 
+enum external_function_type
+{
+  ExternalFunctionWithoutDerivative,
+  ExternalFunctionWithFirstDerivative,
+  ExternalFunctionWithFirstandSecondDerivative,
+  ExternalFunctionNumericalFirstDerivative,
+  ExternalFunctionFirstDerivative,
+  ExternalFunctionNumericalSecondDerivative,
+  ExternalFunctionSecondDerivative
+};
+
 struct Block_contain_type
 {
   int Equation, Variable, Own_Derivative;
@@ -362,6 +383,24 @@ public:
   };
 };
 
+
+class FPUSH_ : public TagWithoutArgument
+{
+public:
+  inline FPUSH_() : TagWithoutArgument(FPUSH)
+  {
+  };
+};
+
+
+class FPOP_ : public TagWithoutArgument
+{
+public:
+  inline FPOP_() : TagWithoutArgument(FPOP)
+  {
+  };
+};
+
 class FDIMT_ : public TagWithOneArgument<unsigned int>
 {
 public:
@@ -734,6 +773,136 @@ public:
   }
 };
 
+
+class FLDTEF_ : public TagWithOneArgument<unsigned int>
+{
+public:
+  inline FLDTEF_() : TagWithOneArgument<unsigned int>::TagWithOneArgument(FLDTEF)
+  {
+  };
+  inline FLDTEF_(unsigned int number) : TagWithOneArgument<unsigned int>::TagWithOneArgument(FLDTEF, number)
+  {
+  };
+  inline unsigned int
+  get_number()
+  {
+    return arg1;
+  }
+};
+
+
+class FSTPTEF_ : public TagWithOneArgument<unsigned int>
+{
+public:
+  inline FSTPTEF_() : TagWithOneArgument<unsigned int>::TagWithOneArgument(FSTPTEF)
+  {
+  };
+  inline FSTPTEF_(unsigned int number) : TagWithOneArgument<unsigned int>::TagWithOneArgument(FSTPTEF, number)
+  {
+  };
+  inline unsigned int
+  get_number()
+  {
+    return arg1;
+  }
+};
+
+class FLDTEFD_ : public TagWithTwoArguments<unsigned int, unsigned int>
+{
+public:
+  inline FLDTEFD_() : TagWithTwoArguments<unsigned int, unsigned int>::TagWithTwoArguments(FLDTEFD)
+  {
+  };
+  inline FLDTEFD_(unsigned int indx, unsigned int row) : TagWithTwoArguments<unsigned int, unsigned int>::TagWithTwoArguments(FLDTEFD, indx, row)
+  {
+  };
+  inline unsigned int
+  get_indx()
+  {
+    return arg1;
+  };
+  inline unsigned int
+  get_row()
+  {
+    return arg2;
+  };
+};
+
+
+class FSTPTEFD_ : public TagWithTwoArguments<unsigned int, unsigned int>
+{
+public:
+  inline FSTPTEFD_() : TagWithTwoArguments<unsigned int, unsigned int>::TagWithTwoArguments(FSTPTEFD)
+  {
+  };
+  inline FSTPTEFD_(unsigned int indx, unsigned int row) : TagWithTwoArguments<unsigned int, unsigned int>::TagWithTwoArguments(FSTPTEFD, indx, row)
+  {
+  };
+  inline unsigned int
+  get_indx()
+  {
+    return arg1;
+  };
+  inline unsigned int
+  get_row()
+  {
+    return arg2;
+  };
+};
+
+class FLDTEFDD_ : public TagWithThreeArguments<unsigned int, unsigned int, unsigned int>
+{
+public:
+  inline FLDTEFDD_() : TagWithThreeArguments<unsigned int, unsigned int, unsigned int>::TagWithThreeArguments(FLDTEFDD)
+  {
+  };
+  inline FLDTEFDD_(unsigned int indx, unsigned int row, unsigned int col) : TagWithThreeArguments<unsigned int, unsigned int, unsigned int>::TagWithThreeArguments(FLDTEFDD, indx, row, col)
+  {
+  };
+  inline unsigned int
+  get_indx()
+  {
+    return arg1;
+  };
+  inline unsigned int
+  get_row()
+  {
+    return arg2;
+  };
+  inline unsigned int
+  get_col()
+  {
+    return arg3;
+  };
+};
+
+
+class FSTPTEFDD_ : public TagWithThreeArguments<unsigned int, unsigned int, unsigned int>
+{
+public:
+  inline FSTPTEFDD_() : TagWithThreeArguments<unsigned int, unsigned int, unsigned int>::TagWithThreeArguments(FSTPTEFDD)
+  {
+  };
+  inline FSTPTEFDD_(unsigned int indx, unsigned int row, unsigned int col) : TagWithThreeArguments<unsigned int, unsigned int, unsigned int>::TagWithThreeArguments(FSTPTEF, indx, row, col)
+  {
+  };
+  inline unsigned int
+  get_indx()
+  {
+    return arg1;
+  };
+  inline unsigned int
+  get_row()
+  {
+    return arg2;
+  };
+  inline unsigned int
+  get_col()
+  {
+    return arg3;
+  };
+};
+
 class FLDVS_ : public TagWithTwoArguments<uint8_t, unsigned int>
 {
 public:
@@ -861,6 +1030,156 @@ public:
   };
 };
 
+
+class FCALL_ : public TagWithFourArguments<unsigned int, unsigned int, string, unsigned int>
+{
+  string func_name;
+  string arg_func_name;
+  unsigned int add_input_arguments, row, col;
+  external_function_type function_type;
+public:
+  inline FCALL_() : TagWithFourArguments<unsigned int, unsigned int, string, unsigned int>::TagWithFourArguments(FCALL)
+  {
+    arg_func_name = "";
+    add_input_arguments = 0;
+    row = 0;
+    col = 0;
+    function_type = ExternalFunctionWithoutDerivative;
+  };
+  inline FCALL_(unsigned int nb_output_arguments, unsigned int nb_input_arguments, string f_name, unsigned int indx) :
+    TagWithFourArguments<unsigned int, unsigned int, string, unsigned int>::TagWithFourArguments(FCALL, nb_output_arguments, nb_input_arguments, f_name, indx)
+  {
+    arg_func_name = "";
+    add_input_arguments = 0;
+    row = 0;
+    col = 0;
+    function_type = ExternalFunctionWithoutDerivative;
+    func_name = f_name;
+  };
+  inline string
+  get_function_name()
+  {
+    //printf("get_function_name => func_name=%s\n",func_name.c_str());fflush(stdout);
+    return func_name;
+  };
+  inline unsigned int
+  get_nb_output_arguments()
+  {
+    return arg1;
+  };
+  inline unsigned int
+  get_nb_input_arguments()
+  {
+    return arg2;
+  };
+  inline unsigned int
+  get_indx()
+  {
+    return arg4;
+  };
+  inline void
+  set_arg_func_name(string arg_arg_func_name)
+  {
+    arg_func_name = arg_arg_func_name;
+  };
+  inline string
+  get_arg_func_name()
+  {
+    return arg_func_name;
+  };
+  inline void
+  set_nb_add_input_arguments(unsigned int arg_add_input_arguments)
+  {
+    add_input_arguments = arg_add_input_arguments;
+  };
+  inline unsigned int
+  get_nb_add_input_arguments()
+  {
+    return add_input_arguments;
+  };
+  inline void
+  set_row(unsigned int arg_row)
+  {
+    row = arg_row;
+  };
+  inline unsigned int
+  get_row()
+  {
+    return row;
+  }
+  inline void
+  set_col(unsigned int arg_col)
+  {
+    col = arg_col;
+  };
+  inline unsigned int
+  get_col()
+  {
+    return col;
+  };
+  inline void
+  set_function_type(external_function_type arg_function_type)
+    {
+      function_type = arg_function_type;
+    };
+  inline external_function_type
+  get_function_type()
+    {
+      return(function_type);
+    }
+  inline void
+  write(ostream &CompileCode, unsigned int &instruction_number)
+  {
+    CompileCode.write(reinterpret_cast<char *>(&op_code), sizeof(op_code));
+    CompileCode.write(reinterpret_cast<char *>(&arg1), sizeof(arg1));
+    CompileCode.write(reinterpret_cast<char *>(&arg2), sizeof(arg2));
+    CompileCode.write(reinterpret_cast<char *>(&arg4), sizeof(arg4));
+    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));
+    int 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 BYTE_CODE
+
+  inline uint8_t *
+  load(uint8_t *code)
+  {
+    op_code = FCALL; code += sizeof(op_code);
+    memcpy(&arg1, code, sizeof(arg1)); code += sizeof(arg1);
+    memcpy(&arg2, code, sizeof(arg2)); code += sizeof(arg2);
+    memcpy(&arg4, code, sizeof(arg4)); code += sizeof(arg4);
+    memcpy(&add_input_arguments, code, sizeof(add_input_arguments)); code += sizeof(add_input_arguments);
+    memcpy(&row, code, sizeof(row)); code += sizeof(row);
+    memcpy(&col, code, sizeof(col)); code += sizeof(col);
+    memcpy(&function_type, code, sizeof(function_type)); code += sizeof(function_type);
+    int size;
+    memcpy(&size, code, sizeof(size)); code += sizeof(size);
+    char* name = (char*)mxMalloc((size+1)*sizeof(char));
+    memcpy(name, code, size); code += size;
+    name[size] = NULL;
+    func_name = name;
+    mxFree(name);
+    memcpy(&size, code, sizeof(size)); code += sizeof(size);
+    name = (char*)mxMalloc((size+1)*sizeof(char));
+    memcpy(name, code, size); code += size;
+    name[size] = NULL;
+    arg_func_name = name;
+    mxFree(name);
+    return code;
+  }
+#endif
+};
+
+
 class FNUMEXPR_ : public TagWithOneArgument<ExpressionType>
 {
 private:
@@ -1473,6 +1792,78 @@ public:
             tags_liste.push_back(make_pair(FJMP, code));
             code += sizeof(FJMP_);
             break;
+          case FCALL:
+            {
+# ifdef DEBUGL
+              mexPrintf("FCALL\n");
+# endif
+              FCALL_ *fcall = new FCALL_;
+
+              code = fcall->load(code);
+
+              tags_liste.push_back(make_pair(FCALL, fcall));
+# ifdef DEBUGL
+              mexPrintf("FCALL finish\n");mexEvalString("drawnow;");
+              mexPrintf("-- *code=%d\n",*code);mexEvalString("drawnow;");
+# endif
+            }
+            break;
+          case FPUSH:
+# ifdef DEBUGL
+            mexPrintf("FPUSH\n");
+# endif
+            tags_liste.push_back(make_pair(FPUSH, code));
+            code += sizeof(FPUSH_);
+            break;
+          case FPOP:
+# ifdef DEBUGL
+            mexPrintf("FPOP\n");
+# endif
+            tags_liste.push_back(make_pair(FPOP, code));
+            code += sizeof(FPOP_);
+            break;
+          case FLDTEF:
+# ifdef DEBUGL
+            mexPrintf("FLDTEF\n");
+# endif
+            tags_liste.push_back(make_pair(FLDTEF, code));
+            code += sizeof(FLDTEF_);
+            break;
+          case FSTPTEF:
+# ifdef DEBUGL
+            mexPrintf("FSTPTEF\n");
+# endif
+            tags_liste.push_back(make_pair(FSTPTEF, code));
+            code += sizeof(FSTPTEF_);
+            break;
+          case FLDTEFD:
+# ifdef DEBUGL
+            mexPrintf("FLDTEFD\n");
+# endif
+            tags_liste.push_back(make_pair(FLDTEFD, code));
+            code += sizeof(FLDTEFD_);
+            break;
+          case FSTPTEFD:
+# ifdef DEBUGL
+            mexPrintf("FSTPTEFD\n");
+# endif
+            tags_liste.push_back(make_pair(FSTPTEFD, code));
+            code += sizeof(FSTPTEFD_);
+            break;
+          case FLDTEFDD:
+# ifdef DEBUGL
+            mexPrintf("FLDTEFDD\n");
+# endif
+            tags_liste.push_back(make_pair(FLDTEFDD, code));
+            code += sizeof(FLDTEFDD_);
+            break;
+          case FSTPTEFDD:
+# ifdef DEBUGL
+            mexPrintf("FSTPTEFDD\n");
+# endif
+            tags_liste.push_back(make_pair(FSTPTEFDD, code));
+            code += sizeof(FSTPTEFDD_);
+            break;
           default:
             mexPrintf("Unknown Tag value=%d code=%x\n", *code, code);
             done = true;
diff --git a/DynamicModel.cc b/DynamicModel.cc
index f071edd866a552c8f9134cceddd86cfcf7acccd7..5013eee791d8ed6f67405a9116a699f64b4fb76a 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -252,6 +252,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
       unsigned int block_size = getBlockSize(block);
       unsigned int block_mfs = getBlockMfs(block);
       unsigned int block_recursive = block_size - block_mfs;
+      deriv_node_temp_terms_t tef_terms;
       /*unsigned int block_exo_size = exo_block[block].size();
       unsigned int block_exo_det_size = exo_det_block[block].size();
       unsigned int block_other_endo_size = other_endo_block[block].size();*/
@@ -475,10 +476,13 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
                    it != v_temporary_terms[block][i].end(); it++)
                 {
+                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
+                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, tef_terms);
+
                   output << "  " <<  sps;
-                  (*it)->writeOutput(output, local_output_type, local_temporary_terms);
+                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, tef_terms);
                   output << " = ";
-                  (*it)->writeOutput(output, local_output_type, tt2);
+                  (*it)->writeOutput(output, local_output_type, tt2, tef_terms);
                   // Insert current node into tt2
                   tt2.insert(*it);
                   output << ";" << endl;
@@ -1070,6 +1074,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
   BinaryOpNode *eq_node;
   Uff Uf[symbol_table.endo_nbr()];
   map<expr_t, int> reference_count;
+  deriv_node_temp_terms_t tef_terms;
   vector<int> feedback_variables;
   bool file_open = false;
 
@@ -1183,9 +1188,12 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
               for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
                    it != v_temporary_terms[block][i].end(); it++)
                 {
+                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
+                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
+
                   FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
                   fnumexpr.write(code_file, instruction_number);
-                  (*it)->compile(code_file, instruction_number, false, tt2, map_idx, true, false);
+                  (*it)->compile(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
                   FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
                   fstpt.write(code_file, instruction_number);
                   // Insert current node into tt2
diff --git a/ExprNode.cc b/ExprNode.cc
index 39356b10ebb75fbdb957922907c0179e3a1addc1..676f0fac3810c8442b83f9a05c61b3347141b49d 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -157,6 +157,16 @@ ExprNode::writeOutput(ostream &output, ExprNodeOutputType output_type, const tem
   writeOutput(output, output_type, temporary_terms, tef_terms);
 }
 
+void
+ExprNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+{
+  deriv_node_temp_terms_t tef_terms;
+  compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
+}
+
+
 void
 ExprNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                       const temporary_terms_t &temporary_terms,
@@ -165,6 +175,16 @@ ExprNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output
   // Nothing to do
 }
 
+
+void
+ExprNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    deriv_node_temp_terms_t &tef_terms) const
+{
+  // Nothing to do
+}
+
 VariableNode *
 ExprNode::createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -294,13 +314,16 @@ NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
 }
 
 double
-NumConstNode::eval(const eval_context_t &eval_context) const throw (EvalException)
+NumConstNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
 {
   return (datatree.num_constants.getDouble(id));
 }
 
 void
-NumConstNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+NumConstNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                      deriv_node_temp_terms_t &tef_terms) const
 {
   FLDC_ fldc(datatree.num_constants.getDouble(id));
   fldc.write(CompileCode, instruction_number);
@@ -710,7 +733,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
 }
 
 double
-VariableNode::eval(const eval_context_t &eval_context) const throw (EvalException)
+VariableNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
 {
   eval_context_t::const_iterator it = eval_context.find(symb_id);
   if (it == eval_context.end())
@@ -720,10 +743,13 @@ VariableNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
 }
 
 void
-VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                      deriv_node_temp_terms_t &tef_terms) const
 {
   if (type == eModelLocalVariable || type == eModFileLocalVariable)
-    datatree.local_variables_table[symb_id]->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+    datatree.local_variables_table[symb_id]->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
   else
     {
       int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
@@ -1593,8 +1619,19 @@ UnaryOpNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType out
   arg->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
 }
 
+void
+UnaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    deriv_node_temp_terms_t &tef_terms) const
+{
+  arg->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+                                     dynamic, steady_dynamic, tef_terms);
+}
+
+
 double
-UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) throw (EvalException)
+UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) throw (EvalException, EvalExternalFunctionException)
 {
   switch (op_code)
     {
@@ -1644,7 +1681,7 @@ UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) throw (EvalException)
 }
 
 double
-UnaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException)
+UnaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
 {
   double v = arg->eval(eval_context);
 
@@ -1652,7 +1689,10 @@ UnaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException
 }
 
 void
-UnaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+UnaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                      deriv_node_temp_terms_t &tef_terms) const
 {
   temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode *>(this));
   if (it != temporary_terms.end())
@@ -1672,10 +1712,10 @@ UnaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, boo
       return;
     }
   if (op_code == oSteadyState)
-    arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, true);
+    arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, true, tef_terms);
   else
     {
-      arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+      arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
       FUNARY_ funary(op_code);
       funary.write(CompileCode, instruction_number);
     }
@@ -2331,7 +2371,7 @@ BinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
 }
 
 double
-BinaryOpNode::eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException)
+BinaryOpNode::eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException, EvalExternalFunctionException)
 {
   switch (op_code)
     {
@@ -2375,7 +2415,7 @@ BinaryOpNode::eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (Eva
 }
 
 double
-BinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException)
+BinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
 {
   double v1 = arg1->eval(eval_context);
   double v2 = arg2->eval(eval_context);
@@ -2384,7 +2424,10 @@ BinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
 }
 
 void
-BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                      deriv_node_temp_terms_t &tef_terms) const
 {
   // If current node is a temporary term
   temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
@@ -2404,8 +2447,8 @@ BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bo
         }
       return;
     }
-  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
-  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
+  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
   FBINARY_ fbinary(op_code);
   fbinary.write(CompileCode, instruction_number);
 }
@@ -2593,6 +2636,18 @@ BinaryOpNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType ou
   arg2->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
 }
 
+void
+BinaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    deriv_node_temp_terms_t &tef_terms) const
+{
+  arg1->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+                                     dynamic, steady_dynamic, tef_terms);
+  arg2->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+                                     dynamic, steady_dynamic, tef_terms);
+}
+
 void
 BinaryOpNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const
 {
@@ -3342,7 +3397,7 @@ TrinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
 }
 
 double
-TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException)
+TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException, EvalExternalFunctionException)
 {
   switch (op_code)
     {
@@ -3356,7 +3411,7 @@ TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v
 }
 
 double
-TrinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException)
+TrinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
 {
   double v1 = arg1->eval(eval_context);
   double v2 = arg2->eval(eval_context);
@@ -3366,7 +3421,10 @@ TrinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalExcepti
 }
 
 void
-TrinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
+TrinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                      bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                      const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                      deriv_node_temp_terms_t &tef_terms) const
 {
   // If current node is a temporary term
   temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<TrinaryOpNode *>(this));
@@ -3386,9 +3444,9 @@ TrinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, b
         }
       return;
     }
-  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
-  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
-  arg3->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
+  arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
+  arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
+  arg3->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic, tef_terms);
   FTRINARY_ ftrinary(op_code);
   ftrinary.write(CompileCode, instruction_number);
 }
@@ -3485,6 +3543,20 @@ TrinaryOpNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType o
   arg3->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
 }
 
+void
+TrinaryOpNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    deriv_node_temp_terms_t &tef_terms) const
+{
+  arg1->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+                                     dynamic, steady_dynamic, tef_terms);
+  arg2->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+                                     dynamic, steady_dynamic, tef_terms);
+  arg3->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+                                     dynamic, steady_dynamic, tef_terms);
+}
+
 void
 TrinaryOpNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const
 {
@@ -3751,8 +3823,12 @@ ExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 expr_t
 ExternalFunctionNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
 {
-  cerr << "ExternalFunctionNode::getChainRuleDerivative: operation impossible!" << endl;
-  exit(EXIT_FAILURE);
+  assert(datatree.external_functions_table.getNargs(symb_id) > 0);
+  vector<expr_t> dargs;
+  for (vector<expr_t>::const_iterator it = arguments.begin();
+     it != arguments.end(); it++)
+    dargs.push_back((*it)->getChainRuleDerivative(deriv_id, recursive_variables));
+  return composeDerivatives(dargs);
 }
 
 void
@@ -3822,6 +3898,56 @@ ExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_typ
   output << "TEF_" << getIndxInTefTerms(symb_id, tef_terms);
 }
 
+unsigned int
+ExternalFunctionNode::compileExternalFunctionArguments(ostream &CompileCode, unsigned int &instruction_number,
+                                                       bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                       const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                       deriv_node_temp_terms_t &tef_terms) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin();
+       it != arguments.end(); it++)
+      (*it)->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx,
+                     dynamic, steady_dynamic, tef_terms);
+  return(arguments.size());
+}
+
+
+void
+ExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                              bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                              const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                              deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<ExternalFunctionNode *>(this));
+  if (it != temporary_terms.end())
+    {
+      if (dynamic)
+        {
+          map_idx_t::const_iterator ii = map_idx.find(idx);
+          FLDT_ fldt(ii->second);
+          fldt.write(CompileCode, instruction_number);
+        }
+      else
+        {
+          map_idx_t::const_iterator ii = map_idx.find(idx);
+          FLDST_ fldst(ii->second);
+          fldst.write(CompileCode, instruction_number);
+        }
+      return;
+    }
+
+  if (!lhs_rhs)
+    {
+      FLDTEF_ fldtef(getIndxInTefTerms(symb_id, tef_terms));
+      fldtef.write(CompileCode, instruction_number);
+    }
+  else
+    {
+      FSTPTEF_ fstptef(getIndxInTefTerms(symb_id, tef_terms));
+      fstptef.write(CompileCode, instruction_number);
+    }
+}
+
 void
 ExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                                   const temporary_terms_t &temporary_terms,
@@ -3899,6 +4025,57 @@ ExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutpu
     }
 }
 
+void
+ExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    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 != eExtFunSetButNoNameProvided);
+
+   for (vector<expr_t>::const_iterator it = arguments.begin();
+       it != arguments.end(); it++)
+    (*it)->compileExternalFunctionOutput(CompileCode, instruction_number, lhs_rhs, temporary_terms,
+                                         map_idx, dynamic, steady_dynamic, tef_terms);
+
+  if (!alreadyWrittenAsTefTerm(symb_id, tef_terms))
+    {
+      tef_terms[make_pair(symb_id, arguments)] = (int) tef_terms.size();
+      int indx = getIndxInTefTerms(symb_id, tef_terms);
+      int second_deriv_symb_id = datatree.external_functions_table.getSecondDerivSymbID(symb_id);
+      assert(second_deriv_symb_id != eExtFunSetButNoNameProvided);
+
+      unsigned int nb_output_arguments = 0;
+      if (symb_id == first_deriv_symb_id &&
+          symb_id == second_deriv_symb_id)
+        nb_output_arguments = 3;
+      else if (symb_id == first_deriv_symb_id)
+        nb_output_arguments = 2;
+      else
+        nb_output_arguments = 1;
+      unsigned int nb_input_arguments = compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs, temporary_terms,
+                                       map_idx, dynamic, steady_dynamic, tef_terms);
+
+      FCALL_ fcall(nb_output_arguments, nb_input_arguments, datatree.symbol_table.getName(symb_id), indx);
+      switch(nb_output_arguments)
+        {
+          case 1:
+            fcall.set_function_type(ExternalFunctionWithoutDerivative);
+            break;
+          case 2:
+            fcall.set_function_type(ExternalFunctionWithFirstDerivative);
+            break;
+          case 3:
+            fcall.set_function_type(ExternalFunctionWithFirstandSecondDerivative);
+            break;
+        }
+      fcall.write(CompileCode, instruction_number);
+      FSTPTEF_ fstptef(indx);
+      fstptef.write(CompileCode, instruction_number);
+    }
+}
+
 void
 ExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                            temporary_terms_t &temporary_terms,
@@ -3907,8 +4084,10 @@ ExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                            vector< vector<temporary_terms_t> > &v_temporary_terms,
                                            int equation) const
 {
-  cerr << "ExternalFunctionNode::computeTemporaryTerms: not implemented" << endl;
-  exit(EXIT_FAILURE);
+  expr_t this2 = const_cast<ExternalFunctionNode *>(this);
+  temporary_terms.insert(this2);
+  first_occurence[this2] = make_pair(Curr_block, equation);
+  v_temporary_terms[Curr_block][equation].insert(this2);
 }
 
 void
@@ -3927,22 +4106,18 @@ ExternalFunctionNode::collectTemporary_terms(const temporary_terms_t &temporary_
     temporary_terms_inuse.insert(idx);
   else
     {
-      //arg->collectTemporary_terms(temporary_terms, result);
+      for (vector<expr_t>::const_iterator it = arguments.begin();
+       it != arguments.end(); it++)
+       (*it)->collectTemporary_terms(temporary_terms, temporary_terms_inuse, Curr_Block);
     }
 }
 
 double
-ExternalFunctionNode::eval(const eval_context_t &eval_context) const throw (EvalException)
+ExternalFunctionNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
 {
-  throw EvalException();
+  throw EvalExternalFunctionException();
 }
 
-void
-ExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
-{
-  cerr << "ExternalFunctionNode::compile: operation impossible!" << endl;
-  exit(EXIT_FAILURE);
-}
 
 pair<int, expr_t>
 ExternalFunctionNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const
@@ -4176,8 +4351,10 @@ FirstDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &referenc
                                                       vector< vector<temporary_terms_t> > &v_temporary_terms,
                                                       int equation) const
 {
-  cerr << "FirstDerivExternalFunctionNode::computeTemporaryTerms: not implemented" << endl;
-  exit(EXIT_FAILURE);
+  expr_t this2 = const_cast<FirstDerivExternalFunctionNode *>(this);
+  temporary_terms.insert(this2);
+  first_occurence[this2] = make_pair(Curr_block, equation);
+  v_temporary_terms[Curr_block][equation].insert(this2);
 }
 
 expr_t
@@ -4233,6 +4410,45 @@ FirstDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
            << LEFT_ARRAY_SUBSCRIPT(output_type) << tmpIndx << RIGHT_ARRAY_SUBSCRIPT(output_type);
 }
 
+void
+FirstDerivExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                              bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                              const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                              deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<FirstDerivExternalFunctionNode *>(this));
+  if (it != temporary_terms.end())
+    {
+      if (dynamic)
+        {
+          map_idx_t::const_iterator ii = map_idx.find(idx);
+          FLDT_ fldt(ii->second);
+          fldt.write(CompileCode, instruction_number);
+        }
+      else
+        {
+          map_idx_t::const_iterator ii = map_idx.find(idx);
+          FLDST_ fldst(ii->second);
+          fldst.write(CompileCode, instruction_number);
+        }
+      return;
+    }
+  int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
+  assert(first_deriv_symb_id != eExtFunSetButNoNameProvided);
+
+  int tmpIndx = inputIndex;
+  if (!lhs_rhs)
+    {
+      FLDTEFD_ fldtefd(getIndxInTefTerms(symb_id, tef_terms), tmpIndx);
+      fldtefd.write(CompileCode, instruction_number);
+    }
+  else
+    {
+      FSTPTEFD_ fstptefd(getIndxInTefTerms(symb_id, tef_terms), tmpIndx);
+      fstptefd.write(CompileCode, instruction_number);
+    }
+}
+
 void
 FirstDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                                             const temporary_terms_t &temporary_terms,
@@ -4326,6 +4542,51 @@ FirstDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Exp
     }
 }
 
+void
+FirstDerivExternalFunctionNode::compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    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 != eExtFunSetButNoNameProvided);
+
+  if (first_deriv_symb_id == symb_id || alreadyWrittenAsTefTerm(first_deriv_symb_id, tef_terms))
+    return;
+
+  unsigned int nb_add_input_arguments = compileExternalFunctionArguments(CompileCode, instruction_number, lhs_rhs, temporary_terms,
+                                       map_idx, dynamic, steady_dynamic, tef_terms);
+  if (first_deriv_symb_id == eExtFunNotSet)
+    {
+      unsigned int nb_input_arguments = 0;
+      unsigned int nb_output_arguments = 1;
+      unsigned int indx = getIndxInTefTerms(symb_id, tef_terms);
+      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(ExternalFunctionNumericalFirstDerivative);
+      fcall.write(CompileCode, instruction_number);
+      FSTPTEFD_ fstptefd(indx, inputIndex);
+      fstptefd.write(CompileCode, instruction_number);
+    }
+  else
+    {
+      tef_terms[make_pair(first_deriv_symb_id, arguments)] = (int) tef_terms.size();
+      int indx = getIndxInTefTerms(symb_id, tef_terms);
+      int second_deriv_symb_id = datatree.external_functions_table.getSecondDerivSymbID(symb_id);
+      assert(second_deriv_symb_id != eExtFunSetButNoNameProvided);
+
+      unsigned int nb_output_arguments = 1;
+
+
+      FCALL_ fcall(nb_output_arguments, nb_add_input_arguments, datatree.symbol_table.getName(first_deriv_symb_id), indx);
+      fcall.set_function_type(ExternalFunctionFirstDerivative);
+      fcall.write(CompileCode, instruction_number);
+    }
+}
+
+
 SecondDerivExternalFunctionNode::SecondDerivExternalFunctionNode(DataTree &datatree_arg,
                                                                  int top_level_symb_id_arg,
                                                                  const vector<expr_t> &arguments_arg,
@@ -4355,8 +4616,10 @@ SecondDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &referen
                                                       vector< vector<temporary_terms_t> > &v_temporary_terms,
                                                       int equation) const
 {
-  cerr << "SecondDerivExternalFunctionNode::computeTemporaryTerms: not implemented" << endl;
-  exit(EXIT_FAILURE);
+  expr_t this2 = const_cast<SecondDerivExternalFunctionNode *>(this);
+  temporary_terms.insert(this2);
+  first_occurence[this2] = make_pair(Curr_block, equation);
+  v_temporary_terms[Curr_block][equation].insert(this2);
 }
 
 expr_t
diff --git a/ExprNode.hh b/ExprNode.hh
index de71231eacb932f11ddc0ded7cdc5726ce961b3c..9aa0d84b3b1cd60f4107c5a7dc44b68f8c9014fe 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -196,6 +196,11 @@ public:
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
 
+  virtual void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    deriv_node_temp_terms_t &tef_terms) const;
+
   //! Computes the set of all variables of a given symbol type in the expression
   /*!
     Variables are stored as integer pairs of the form (symb_id, lag).
@@ -241,8 +246,14 @@ public:
   {
   };
 
-  virtual double eval(const eval_context_t &eval_context) const throw (EvalException) = 0;
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const = 0;
+  class EvalExternalFunctionException : public EvalException
+
+  {
+  };
+
+  virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException ) = 0;
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const = 0;
+  void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
   //! Creates a static version of this node
   /*!
     This method duplicates the current node by creating a similar node from which all leads/lags have been stripped,
@@ -403,8 +414,8 @@ public:
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
-  virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+  virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException );
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
   virtual expr_t toStatic(DataTree &static_datatree) const;
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
@@ -449,8 +460,8 @@ public:
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
-  virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+  virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException );
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
   virtual expr_t toStatic(DataTree &static_datatree) const;
   SymbolType
   get_type() const
@@ -506,6 +517,10 @@ public:
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
+  virtual void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                             bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                             const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                             deriv_node_temp_terms_t &tef_terms) const;
   virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      map<expr_t, pair<int, int> > &first_occurence,
@@ -514,9 +529,9 @@ public:
                                      int equation) const;
   virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
-  static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException);
-  virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+  static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException, EvalExternalFunctionException );
+  virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException );
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
   //! Returns operand
   expr_t
   get_arg() const
@@ -573,6 +588,10 @@ public:
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
+  virtual void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                             bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                             const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                             deriv_node_temp_terms_t &tef_terms) const;
   virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      map<expr_t, pair<int, int> > &first_occurence,
@@ -581,9 +600,9 @@ public:
                                      int equation) const;
   virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
-  static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
-  virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+  static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException, EvalExternalFunctionException );
+  virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException );
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
   virtual expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
   //! Returns first operand
   expr_t
@@ -648,6 +667,10 @@ public:
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
+  virtual void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                             bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                             const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                             deriv_node_temp_terms_t &tef_terms) const;
   virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      map<expr_t, pair<int, int> > &first_occurence,
@@ -656,9 +679,9 @@ public:
                                      int equation) const;
   virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
-  static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
-  virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+  static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException, EvalExternalFunctionException );
+  virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException );
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
   virtual expr_t toStatic(DataTree &static_datatree) const;
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
@@ -711,6 +734,10 @@ public:
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
+  virtual void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    deriv_node_temp_terms_t &tef_terms) const;
   virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      map<expr_t, pair<int, int> > &first_occurence,
@@ -719,8 +746,13 @@ public:
                                      int equation) const;
   virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
-  virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
-  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
+  virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException );
+  unsigned int compileExternalFunctionArguments(ostream &CompileCode, unsigned int &instruction_number,
+                                                       bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                       const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                       deriv_node_temp_terms_t &tef_terms) const;
+
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic, deriv_node_temp_terms_t &tef_terms) const;
   virtual expr_t toStatic(DataTree &static_datatree) const;
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
@@ -763,9 +795,17 @@ public:
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
+                              bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                              const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                              deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
+  virtual void compileExternalFunctionOutput(ostream &CompileCode, unsigned int &instruction_number,
+                                                    bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                                                    const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                                                    deriv_node_temp_terms_t &tef_terms) const;
 };
 
 class SecondDerivExternalFunctionNode : public ExternalFunctionNode
diff --git a/ModFile.cc b/ModFile.cc
index 30d813a4197aa9effe5abc3753e5531ea3d8f059..3c603959aeccb51eb4edb608a3dadc604f798f05 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -147,12 +147,6 @@ ModFile::checkPass()
       exit(EXIT_FAILURE);
     }
 
-  if (byte_code && (external_functions_table.get_total_number_of_unique_model_block_external_functions() > 0))
-    {
-      cerr << "ERROR: In 'model' block, use of external functions is not compatible with 'bytecode'" << endl;
-      exit(EXIT_FAILURE);
-    }
-
   if (mod_file_struct.dsge_var_estimated)
     if (!mod_file_struct.dsge_prior_weight_in_estimated_params)
       {
diff --git a/ModelTree.cc b/ModelTree.cc
index 5ce9b8153205c2f121a3cb9c9e0fbe9757543195..1d158f818f4f1ecca2315886b6c1a1d3f8c7b886 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -264,6 +264,10 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
             {
               val = Id->eval(eval_context);
             }
+          catch (ExprNode::EvalExternalFunctionException &e)
+            {
+              val = 1;
+            }
           catch (ExprNode::EvalException &e)
             {
               cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
@@ -271,7 +275,6 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
               cerr << endl;
               exit(EXIT_FAILURE);
             }
-
           if (fabs(val) < cutoff)
             {
               if (verbose)
@@ -1115,12 +1118,19 @@ ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_n
 {
   // Local var used to keep track of temp nodes already written
   temporary_terms_t tt2;
+  // To store the functions that have already been written in the form TEF* = ext_fun();
+  deriv_node_temp_terms_t tef_terms;
   for (temporary_terms_t::const_iterator it = tt.begin();
        it != tt.end(); it++)
     {
+      if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
+        {
+          (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
+        }
+
       FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
       fnumexpr.write(code_file, instruction_number);
-      (*it)->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic);
+      (*it)->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
       if (dynamic)
         {
           FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index 012fa72e64ad774431c3fb986da2d91c0b9e3a1d..ff372354d4d86cac8620d9e37e6274e18341201c 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -1865,16 +1865,41 @@ ParsingDriver::add_model_var_or_external_function(string *function_name, bool in
               double model_var_arg_dbl;
               if (unaryNode == NULL)
                 {
-                  model_var_arg = (int)numNode->eval(ectmp);
-                  model_var_arg_dbl = numNode->eval(ectmp);
+                  try
+                    {
+                      model_var_arg = (int)numNode->eval(ectmp);
+                    }
+                  catch (ExprNode::EvalException &e)
+                    {
+                    }
+                  try
+                    {
+                      model_var_arg_dbl = numNode->eval(ectmp);
+                    }
+                  catch (ExprNode::EvalException &e)
+                    {
+                    }
+
                 }
               else
                 if (unaryNode->get_op_code() != oUminus)
                   error("A model variable is being treated as if it were a function (i.e., takes an argument that is not an integer).");
                 else
                   {
-                    model_var_arg = (int)unaryNode->eval(ectmp);
-                    model_var_arg_dbl = unaryNode->eval(ectmp);
+                    try
+                      {
+                        model_var_arg = (int)unaryNode->eval(ectmp);
+                      }
+                    catch (ExprNode::EvalException &e)
+                      {
+                      }
+                    try
+                      {
+                        model_var_arg_dbl = unaryNode->eval(ectmp);
+                      }
+                    catch (ExprNode::EvalException &e)
+                      {
+                      }
                   }
 
               if ((double) model_var_arg != model_var_arg_dbl) //make 100% sure int cast didn't lose info
diff --git a/StaticModel.cc b/StaticModel.cc
index 46a536427bb8cbee5018d44eec49c43fd759ef34..9488efda7d828778eed904c5b1ba157093bd388c 100644
--- a/StaticModel.cc
+++ b/StaticModel.cc
@@ -133,14 +133,9 @@ StaticModel::computeTemporaryTermsOrdered()
           expr_t id = it->second.second;
           id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  block_size-1);
         }
-      //temporary_terms_g.insert(temporary_terms_l.begin(), temporary_terms_l.end());
       set<int> temporary_terms_in_use;
       temporary_terms_in_use.clear();
       v_temporary_terms_inuse[block] = temporary_terms_in_use;
-      /*for (int i = 0; i < (int) block_size; i++)
-        for (temporary_terms_t::const_iterator it = v_temporary_terms_local[block][i].begin();
-             it != v_temporary_terms_local[block][i].end(); it++)
-          (*it)->collectTemporary_terms(temporary_terms_g, temporary_terms_in_use, block);*/
       computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]);
     }
 
@@ -222,6 +217,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
   ofstream  output;
   int nze;
   vector<int> feedback_variables;
+  deriv_node_temp_terms_t tef_terms;
   ExprNodeOutputType local_output_type;
 
   local_output_type = oMatlabStaticModelSparse;
@@ -304,10 +300,13 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
               for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
                    it != v_temporary_terms[block][i].end(); it++)
                 {
+                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
+                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, tef_terms);
+
                   output << "  " <<  sps;
-                  (*it)->writeOutput(output, local_output_type, local_temporary_terms);
+                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, tef_terms);
                   output << " = ";
-                  (*it)->writeOutput(output, local_output_type, tt2);
+                  (*it)->writeOutput(output, local_output_type, tt2, tef_terms);
                   // Insert current node into tt2
                   tt2.insert(*it);
                   output << ";" << endl;
@@ -602,6 +601,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
   Uff Uf[symbol_table.endo_nbr()];
   map<expr_t, int> reference_count;
   vector<int> feedback_variables;
+  deriv_node_temp_terms_t tef_terms;
   bool file_open = false;
 
   string main_name = file_name;
@@ -655,11 +655,6 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
 
       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;
       for (i = 0; i < (int) block_size; i++)
         {
           //The Temporary terms
@@ -670,9 +665,12 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
               for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
                    it != v_temporary_terms[block][i].end(); it++)
                 {
+                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
+                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
+
                   FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
                   fnumexpr.write(code_file, instruction_number);
-                  (*it)->compile(code_file, instruction_number, false, tt2, map_idx, false, false);
+                  (*it)->compile(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
                   FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second));
                   fstpst.write(code_file, instruction_number);
                   // Insert current node into tt2
@@ -738,6 +736,13 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
         }
       FENDEQU_ fendequ;
       fendequ.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;
+
       // The Jacobian if we have to solve the block
       if    (simulation_type != EVALUATE_BACKWARD
              && simulation_type != EVALUATE_FORWARD)
@@ -858,9 +863,14 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
               for (temporary_terms_t::const_iterator it = v_temporary_terms_local[block][i].begin();
                    it != v_temporary_terms_local[block][i].end(); it++)
                 {
+                  if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
+                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);
+
                   FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx2[block].find((*it)->idx)->second));
                   fnumexpr.write(code_file, instruction_number);
-                  (*it)->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false);
+
+                  (*it)->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);
+
                   FSTPST_ fstpst((int)(map_idx2[block].find((*it)->idx)->second));
                   fstpst.write(code_file, instruction_number);
                   // Insert current node into tt2
@@ -937,7 +947,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
             FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
             fnumexpr.write(code_file, instruction_number);
           }
-          compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx2[block], tt2);
+          compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx2[block], tt2/*temporary_terms*/);
           {
             FSTPG2_ fstpg2(0,0);
             fstpg2.write(code_file, instruction_number);
@@ -957,7 +967,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
               FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0);
               fnumexpr.write(code_file, instruction_number);
 
-              compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx2[block], tt2);
+              compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx2[block], tt2/*temporary_terms*/);
 
               FSTPG2_ fstpg2(eq,var);
               fstpg2.write(code_file, instruction_number);