diff --git a/DynamicModel.cc b/DynamicModel.cc
index 2d53e19aa00b9ed03ecb192be8ca671337ec537d..8ec2e662feed459da5c517f95a53b1a373d61803 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -2123,11 +2123,13 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
                                     julia ? oJuliaDynamicModel : oMatlabDynamicModel);
 
   deriv_node_temp_terms_t tef_terms;
+  temporary_terms_t temp_term_empty;
   temporary_terms_t temp_term_union = temporary_terms_res;
+  temporary_terms_t temp_term_union_m_1;
 
   writeModelLocalVariables(model_local_vars_output, output_type, tef_terms);
 
-  writeTemporaryTerms(temporary_terms_res, model_output, output_type, tef_terms);
+  writeTemporaryTerms(temporary_terms_res, temp_term_union_m_1, model_output, output_type, tef_terms);
 
   writeModelEquations(model_output, output_type);
 
@@ -2135,12 +2137,13 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
   int hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr;
 
   // Writing Jacobian
+  temp_term_union_m_1 = temp_term_union;
   temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
   if (!first_derivatives.empty())
     if (julia)
-      writeTemporaryTerms(temp_term_union, jacobian_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_empty, jacobian_output, output_type, tef_terms);
     else
-      writeTemporaryTerms(temporary_terms_g1, jacobian_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_union_m_1, jacobian_output, output_type, tef_terms);
   for (first_derivatives_t::const_iterator it = first_derivatives.begin();
        it != first_derivatives.end(); it++)
     {
@@ -2155,12 +2158,13 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
     }
 
   // Writing Hessian
+  temp_term_union_m_1 = temp_term_union;
   temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
   if (!second_derivatives.empty())
     if (julia)
-      writeTemporaryTerms(temp_term_union, hessian_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_empty, hessian_output, output_type, tef_terms);
     else
-      writeTemporaryTerms(temporary_terms_g2, hessian_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_union_m_1, hessian_output, output_type, tef_terms);
   int k = 0; // Keep the line of a 2nd derivative in v2
   for (second_derivatives_t::const_iterator it = second_derivatives.begin();
        it != second_derivatives.end(); it++)
@@ -2223,12 +2227,13 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
     }
 
   // Writing third derivatives
+  temp_term_union_m_1 = temp_term_union;
   temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
   if (!third_derivatives.empty())
     if (julia)
-      writeTemporaryTerms(temp_term_union, third_derivatives_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_empty, third_derivatives_output, output_type, tef_terms);
     else
-      writeTemporaryTerms(temporary_terms_g3, third_derivatives_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_union_m_1, third_derivatives_output, output_type, tef_terms);
   k = 0; // Keep the line of a 3rd derivative in v3
   for (third_derivatives_t::const_iterator it = third_derivatives.begin();
        it != third_derivatives.end(); it++)
@@ -4036,7 +4041,8 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
   deriv_node_temp_terms_t tef_terms;
   writeModelLocalVariables(paramsDerivsFile, output_type, tef_terms);
 
-  writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, output_type, tef_terms);
+  temporary_terms_t temp_terms_empty;
+  writeTemporaryTerms(params_derivs_temporary_terms, temp_terms_empty, paramsDerivsFile, output_type, tef_terms);
 
   // Write parameter derivative
   paramsDerivsFile << "rp = zeros(" << equation_number() << ", "
diff --git a/ModelTree.cc b/ModelTree.cc
index 1fa26332bd651f013cc390fc242cd86e294a9c42..b7fabb537d529f2ae6de4efec1f3f1aeeba6ac05 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -1262,33 +1262,34 @@ ModelTree::computeTemporaryTerms(bool is_matlab)
 }
 
 void
-ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, ostream &output,
+ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, const temporary_terms_t &ttm1, ostream &output,
                                ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const
 {
   // Local var used to keep track of temp nodes already written
-  temporary_terms_t tt2;
+  temporary_terms_t tt2 = ttm1;
   for (temporary_terms_t::const_iterator it = tt.begin();
        it != tt.end(); it++)
-    {
-      if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
-        (*it)->writeExternalFunctionOutput(output, output_type, tt2, tef_terms);
+    if (ttm1.find(*it) == ttm1.end())
+      {
+        if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
+          (*it)->writeExternalFunctionOutput(output, output_type, tt2, tef_terms);
 
-      if (IS_C(output_type))
-        output << "double ";
-      else if (IS_JULIA(output_type))
-        output << "  @inbounds const ";
+        if (IS_C(output_type))
+          output << "double ";
+        else if (IS_JULIA(output_type))
+          output << "  @inbounds const ";
 
-      (*it)->writeOutput(output, output_type, tt, tef_terms);
-      output << " = ";
-      (*it)->writeOutput(output, output_type, tt2, tef_terms);
+        (*it)->writeOutput(output, output_type, tt, tef_terms);
+        output << " = ";
+        (*it)->writeOutput(output, output_type, tt2, tef_terms);
 
-      if (IS_C(output_type) || IS_MATLAB(output_type))
-        output << ";";
-      output << endl;
+        if (IS_C(output_type) || IS_MATLAB(output_type))
+          output << ";";
+        output << endl;
 
-      // Insert current node into tt2
-      tt2.insert(*it);
-    }
+        // Insert current node into tt2
+        tt2.insert(*it);
+      }
 }
 
 void
diff --git a/ModelTree.hh b/ModelTree.hh
index 4db8fb8673999f8ecf1afbb7df371d9f5958067e..ee4e73cf64bbad6ca26f81ce14e84334ae141ea7 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -185,7 +185,7 @@ protected:
   //! Computes temporary terms for the file containing parameters derivatives
   void computeParamsDerivativesTemporaryTerms();
 //! Writes temporary terms
-  void writeTemporaryTerms(const temporary_terms_t &tt, ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
+  void writeTemporaryTerms(const temporary_terms_t &tt, const temporary_terms_t &ttm1, ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
   //! Compiles temporary terms
   void compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const;
   //! Adds informations for simulation in a binary file
diff --git a/StaticModel.cc b/StaticModel.cc
index d43c09a85aa8efef288570f08b4ec853ba13b70d..e2eb97d0a9a898efe4f6f9f6f0c588fa4a356b2c 100644
--- a/StaticModel.cc
+++ b/StaticModel.cc
@@ -1212,11 +1212,13 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
                                     julia ? oJuliaStaticModel : oMatlabStaticModel);
 
   deriv_node_temp_terms_t tef_terms;
+  temporary_terms_t temp_term_empty;
   temporary_terms_t temp_term_union = temporary_terms_res;
+  temporary_terms_t temp_term_union_m_1;
 
   writeModelLocalVariables(model_local_vars_output, output_type, tef_terms);
 
-  writeTemporaryTerms(temporary_terms_res, model_output, output_type, tef_terms);
+  writeTemporaryTerms(temporary_terms_res, temp_term_union_m_1, model_output, output_type, tef_terms);
 
   writeModelEquations(model_output, output_type);
 
@@ -1225,12 +1227,13 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
   int hessianColsNbr = JacobianColsNbr*JacobianColsNbr;
 
   // Write Jacobian w.r. to endogenous only
+  temp_term_union_m_1 = temp_term_union;
   temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
   if (!first_derivatives.empty())
     if (julia)
-      writeTemporaryTerms(temp_term_union, jacobian_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_empty, jacobian_output, output_type, tef_terms);
     else
-      writeTemporaryTerms(temporary_terms_g1, jacobian_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_union_m_1, jacobian_output, output_type, tef_terms);
   for (first_derivatives_t::const_iterator it = first_derivatives.begin();
        it != first_derivatives.end(); it++)
     {
@@ -1246,12 +1249,13 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
 
   int g2ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
   // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
+  temp_term_union_m_1 = temp_term_union;
   temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
   if (!second_derivatives.empty())
     if (julia)
-      writeTemporaryTerms(temp_term_union, hessian_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_empty, hessian_output, output_type, tef_terms);
     else
-      writeTemporaryTerms(temporary_terms_g2, hessian_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_union_m_1, hessian_output, output_type, tef_terms);
   int k = 0; // Keep the line of a 2nd derivative in v2
   for (second_derivatives_t::const_iterator it = second_derivatives.begin();
        it != second_derivatives.end(); it++)
@@ -1313,12 +1317,13 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
     }
 
   // Writing third derivatives
+  temp_term_union_m_1 = temp_term_union;
   temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
   if (!third_derivatives.empty())
     if (julia)
-      writeTemporaryTerms(temp_term_union, third_derivatives_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_empty, third_derivatives_output, output_type, tef_terms);
     else
-      writeTemporaryTerms(temporary_terms_g3, third_derivatives_output, output_type, tef_terms);
+      writeTemporaryTerms(temp_term_union, temp_term_union_m_1, third_derivatives_output, output_type, tef_terms);
   k = 0; // Keep the line of a 3rd derivative in v3
   for (third_derivatives_t::const_iterator it = third_derivatives.begin();
        it != third_derivatives.end(); it++)
@@ -2207,7 +2212,8 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
   deriv_node_temp_terms_t tef_terms;
   writeModelLocalVariables(paramsDerivsFile, output_type, tef_terms);
 
-  writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, output_type, tef_terms);
+  temporary_terms_t temp_terms_empty;
+  writeTemporaryTerms(params_derivs_temporary_terms, temp_terms_empty, paramsDerivsFile, output_type, tef_terms);
 
   // Write parameter derivative
   paramsDerivsFile << "rp = zeros(" << equation_number() << ", "