diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 7bcf7a6e10a36a93a8ea59ef7bb8058e592acede..5a6fb32c08f9762b6c736b3c0dd21b260d3a6a81 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -449,6 +449,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
         else
           sps = "";
       // The equations
+      temporary_terms_idxs_t temporary_terms_idxs;
       for (unsigned int i = 0; i < block_size; i++)
         {
           temporary_terms_t tt2;
@@ -460,12 +461,12 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
                    it != v_temporary_terms[block][i].end(); it++)
                 {
                   if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
-                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, tef_terms);
+                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
 
                   output << "  " <<  sps;
-                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, tef_terms);
+                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, temporary_terms_idxs, tef_terms);
                   output << " = ";
-                  (*it)->writeOutput(output, local_output_type, tt2, tef_terms);
+                  (*it)->writeOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
                   // Insert current node into tt2
                   tt2.insert(*it);
                   output << ";" << endl;
@@ -1522,49 +1523,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
 void
 DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
 {
-  string filename = dynamic_basename + ".m";
-
-  ofstream mDynamicModelFile;
-  mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
-  if (!mDynamicModelFile.is_open())
-    {
-      cerr << "Error: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-  mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x, params, steady_state, it_)" << endl
-                    << "%" << endl
-                    << "% Status : Computes dynamic model for Dynare" << endl
-                    << "%" << endl
-                    << "% Inputs :" << endl
-                    << "%   y         [#dynamic variables by 1] double    vector of endogenous variables in the order stored" << endl
-                    << "%                                                 in M_.lead_lag_incidence; see the Manual" << endl
-                    << "%   x         [nperiods by M_.exo_nbr] double     matrix of exogenous variables (in declaration order)" << endl
-                    << "%                                                 for all simulation periods" << endl
-                    << "%   steady_state  [M_.endo_nbr by 1] double       vector of steady state values" << endl
-                    << "%   params    [M_.param_nbr by 1] double          vector of parameter values in declaration order" << endl
-                    << "%   it_       scalar double                       time period for exogenous variables for which to evaluate the model" << endl
-                    << "%" << endl
-                    << "% Outputs:" << endl
-                    << "%   residual  [M_.endo_nbr by 1] double    vector of residuals of the dynamic model equations in order of " << endl
-                    << "%                                          declaration of the equations." << endl
-                    << "%                                          Dynare may prepend auxiliary equations, see M_.aux_vars" << endl
-                    << "%   g1        [M_.endo_nbr by #dynamic variables] double    Jacobian matrix of the dynamic model equations;" << endl
-                    << "%                                                           rows: equations in order of declaration" << endl
-                    << "%                                                           columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names" << endl
-                    << "%   g2        [M_.endo_nbr by (#dynamic variables)^2] double   Hessian matrix of the dynamic model equations;" << endl
-                    << "%                                                              rows: equations in order of declaration" << endl
-                    << "%                                                              columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names" << endl
-                    << "%   g3        [M_.endo_nbr by (#dynamic variables)^3] double   Third order derivative matrix of the dynamic model equations;" << endl
-                    << "%                                                              rows: equations in order of declaration" << endl
-                    << "%                                                              columns: variables in order stored in M_.lead_lag_incidence followed by the ones in M_.exo_names" << endl
-                    << "%" << endl
-                    << "%" << endl
-                    << "% Warning : this file is generated automatically by Dynare" << endl
-                    << "%           from model file (.mod)" << endl << endl;
-
-  writeDynamicModel(mDynamicModelFile, false, false);
-  mDynamicModelFile << "end" << endl; // Close *_dynamic function
-  mDynamicModelFile.close();
+  writeDynamicModel(dynamic_basename, false, false);
 }
 
 void
@@ -1600,26 +1559,7 @@ DynamicModel::writeVarExpectationCalls(ostream &output) const
 void
 DynamicModel::writeDynamicJuliaFile(const string &basename) const
 {
-  string filename = basename + "Dynamic.jl";
-
-  ofstream output;
-  output.open(filename.c_str(), ios::out | ios::binary);
-  if (!output.is_open())
-    {
-      cerr << "Error: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  output << "module " << basename << "Dynamic" << endl
-         << "#" << endl
-         << "# NB: this file was automatically generated by Dynare" << endl
-         << "#     from " << basename << ".mod" << endl
-         << "#" << endl
-         << "using Utils" << endl << endl
-         << "export dynamic!" << endl << endl;
-  writeDynamicModel(output, false, true);
-  output << "end" << endl;
-  output.close();
+  writeDynamicModel(basename, false, true);
 }
 
 void
@@ -2201,200 +2141,374 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
   chdir("..");
 }
 
+void
+DynamicModel::writeWrapperFunctions(const string &basename, const string &ending) const
+{
+  string name;
+  if (ending == "g1")
+    name = basename + "_resid_g1";
+  else if (ending == "g2")
+    name= basename + "_resid_g1_g2";
+  else if (ending == "g3")
+    name = basename + "_resid_g1_g2_g3";
+
+  string filename = name + ".m";
+  ofstream output;
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  if (ending == "g1")
+    output << "function [residual, g1] = " << name << "(T, y, x, params, steady_state, it_)" << endl
+           << "% function [residual, g1] = " << name << "(T, y, x, params, steady_state, it_)" << endl;
+  else if (ending == "g2")
+    output << "function [residual, g1, g2] = " << name << "(T, y, x, params, steady_state, it_)" << endl
+           << "% function [residual, g1, g2] = " << name << "(T, y, x, params, steady_state, it_)" << endl;
+  else if (ending == "g3")
+    output << "function [residual, g1, g2, g3] = " << name << "(T, y, x, params, steady_state, it_)" << endl
+           << "% function [residual, g1, g2, g3] = " << name << "(T, y, x, params, steady_state, it_)" << endl;
+
+  output << "%" << endl
+         << "% Wrapper function automatically created by Dynare" << endl
+         << "%" << endl << endl;
+
+  if (ending == "g1")
+    output << "    T = " << basename + "_" + ending + "_tt(T, y, x, params, steady_state, it_);" << endl
+           << "    residual = " << basename + "_resid(T, y, x, params, steady_state, it_, false);" << endl
+           << "    g1       = " << basename + "_g1(T, y, x, params, steady_state, it_, false);" << endl;
+  else if (ending == "g2")
+    output << "    T = " << basename + "_" + ending + "_tt(T, y, x, params, steady_state, it_);" << endl
+           << "    [residual, g1] = " << basename + "_resid_g1(T, y, x, params, steady_state, it_, false);" << endl
+           << "    g2       = " << basename + "_g2(T, y, x, params, steady_state, it_, false);" << endl;
+  else if (ending == "g3")
+    output << "    T = " << basename + "_" + ending + "_tt(T, y, x, params, steady_state, it_);" << endl
+           << "    [residual, g1, g2] = " << basename + "_resid_g1(T, y, x, params, steady_state, it_, false);" << endl
+           << "    g3       = " << basename + "_g3(T, y, x, params, steady_state, it_, false);" << endl;
+
+  output << endl << "end" << endl;
+  output.close();
+}
+
+void
+DynamicModel::writeDynamicModelHelper(const string &name, const string &retvalname,
+                                      const string &name_tt, size_t ttlen,
+                                      const string &previous_tt_name,
+                                      const ostringstream &init_s,
+                                      const ostringstream &end_s,
+                                      const ostringstream &s, const ostringstream &s_tt) const
+{
+  string filename =  name_tt + ".m";
+  ofstream output;
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  output << "function T = " << name_tt << "(T, y, x, params, steady_state, it_)" << endl
+         << "% function T = " << name_tt << "(T, y, x, params, steady_state, it_)" << endl
+         << "%" << endl
+         << "% File created by Dynare Preprocessor from .mod file" << endl
+         << "%" << endl
+         << "% Inputs:" << endl
+         << "%   T             [#temp variables by 1]     double  vector of temporary terms to be filled by function" << endl
+         << "%   y             [#dynamic variables by 1]  double  vector of endogenous variables in the order stored" << endl
+         << "%                                                    in M_.lead_lag_incidence; see the Manual" << endl
+         << "%   x             [nperiods by M_.exo_nbr]   double  matrix of exogenous variables (in declaration order)" << endl
+         << "%                                                    for all simulation periods" << endl
+         << "%   steady_state  [M_.endo_nbr by 1]         double  vector of steady state values" << endl
+         << "%   params        [M_.param_nbr by 1]        double  vector of parameter values in declaration order" << endl
+         << "%   it_           scalar                     double  time period for exogenous variables for which" << endl
+         << "%                                                    to evaluate the model" << endl
+         << "%" << endl
+         << "% Output:" << endl
+         << "%   T           [#temp variables by 1]       double  vector of temporary terms" << endl
+         << "%" << endl << endl
+         << "assert(length(T) >= " << ttlen << ");" << endl
+         << endl;
+
+  if (!previous_tt_name.empty())
+    output << "T = " << previous_tt_name << "(T, y, x, params, steady_state, it_);" << endl << endl;
+
+  output << s_tt.str() << endl
+         << "end" << endl;
+  output.close();
+
+  filename = name + ".m";
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  output << "function " << retvalname << " = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
+         << "% function " << retvalname << " = " << name << "(T, y, x, params, steady_state, it_, T_flag)" << endl
+         << "%" << endl
+         << "% File created by Dynare Preprocessor from .mod file" << endl
+         << "%" << endl
+         << "% Inputs:" << endl
+         << "%   T             [#temp variables by 1]     double   vector of temporary terms to be filled by function" << endl
+         << "%   y             [#dynamic variables by 1]  double   vector of endogenous variables in the order stored" << endl
+         << "%                                                     in M_.lead_lag_incidence; see the Manual" << endl
+         << "%   x             [nperiods by M_.exo_nbr]   double   matrix of exogenous variables (in declaration order)" << endl
+         << "%                                                     for all simulation periods" << endl
+         << "%   steady_state  [M_.endo_nbr by 1]         double   vector of steady state values" << endl
+         << "%   params        [M_.param_nbr by 1]        double   vector of parameter values in declaration order" << endl
+         << "%   it_           scalar                     double   time period for exogenous variables for which" << endl
+         << "%                                                     to evaluate the model" << endl
+         << "%   T_flag        boolean                    boolean  flag saying whether or not to calculate temporary terms" << endl
+         << "%" << endl
+         << "% Output:" << endl
+         << "%   " << retvalname << endl
+         << "%" << endl << endl;
+
+  if (!name_tt.empty())
+    output << "if T_flag" << endl
+           << "    T = " << name_tt << "(T, y, x, params, steady_state, it_);" << endl
+           << "end" << endl;
+
+  output << init_s.str() << endl
+         << s.str()
+         << end_s.str() << endl
+         << "end" << endl;
+  output.close();
+}
+
 void
 DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const
 {
-  ostringstream model_local_vars_output;  // Used for storing model local vars
-  ostringstream model_output;             // Used for storing model temp vars and equations
-  ostringstream jacobian_output;          // Used for storing jacobian equations
-  ostringstream hessian_output;           // Used for storing Hessian equations
-  ostringstream third_derivatives_output; // Used for storing third order derivatives equations
+  writeDynamicModel("", DynamicOutput, use_dll, julia);
+}
+
+void
+DynamicModel::writeDynamicModel(const string &dynamic_basename, bool use_dll, bool julia) const
+{
+  ofstream DynamicOutput;
+  writeDynamicModel(dynamic_basename, DynamicOutput, use_dll, julia);
+}
+
+void
+DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &DynamicOutput, bool use_dll, bool julia) const
+{
+  ostringstream model_tt_output;             // Used for storing model temp vars
+  ostringstream model_output;                // Used for storing model equations
+  ostringstream jacobian_tt_output;          // Used for storing jacobian temp vars
+  ostringstream jacobian_output;             // Used for storing jacobian equations
+  ostringstream hessian_tt_output;           // Used for storing Hessian temp vars
+  ostringstream hessian_output;              // Used for storing Hessian equations
+  ostringstream third_derivatives_tt_output; // Used for storing third order derivatives temp terms
+  ostringstream third_derivatives_output;    // Used for storing third order derivatives equations
 
   ExprNodeOutputType output_type = (use_dll ? oCDynamicModel :
                                     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, temp_term_union_m_1, model_output, output_type, tef_terms);
-
-  writeModelEquations(model_output, output_type);
+  temporary_terms_t temp_term_union;
+  temporary_terms_idxs_t temp_term_union_idxs;
+
+  for (map<expr_t, expr_t>::const_iterator it = temporary_terms_mlv.begin();
+       it != temporary_terms_mlv.end(); it++)
+    temp_term_union.insert(it->first);
+  writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_mlv,
+                                        temporary_terms_mlv_idxs,
+                                        model_tt_output, output_type, tef_terms);
+  temp_term_union_idxs = temporary_terms_mlv_idxs;
+
+  writeTemporaryTerms(temporary_terms_res, temporary_terms_res_idxs,
+                      temp_term_union, temp_term_union_idxs,
+                      model_tt_output, output_type, tef_terms);
+  temp_term_union.insert(temporary_terms_res.begin(), temporary_terms_res.end());
+  temp_term_union_idxs.insert(temp_term_union_idxs.end(),
+                              temporary_terms_res_idxs.begin(), temporary_terms_res_idxs.end());
+
+  writeModelEquations(model_output, output_type, temp_term_union, temp_term_union_idxs);
 
   int nrows = equations.size();
   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, temp_term_empty, jacobian_output, output_type, tef_terms);
-    else
-      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++)
     {
-      int eq = it->first.first;
-      int var = it->first.second;
-      expr_t d1 = it->second;
+      writeTemporaryTerms(temporary_terms_g1, temporary_terms_g1_idxs,
+                          temp_term_union, temp_term_union_idxs,
+                          jacobian_tt_output, output_type, tef_terms);
+      temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
+      temp_term_union_idxs.insert(temp_term_union_idxs.end(),
+                                  temporary_terms_g1_idxs.begin(), temporary_terms_g1_idxs.end());
 
-      jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
-      jacobian_output << "=";
-      d1->writeOutput(jacobian_output, output_type, temp_term_union, tef_terms);
-      jacobian_output << ";" << endl;
+      for (first_derivatives_t::const_iterator it = first_derivatives.begin();
+           it != first_derivatives.end(); it++)
+        {
+          int eq = it->first.first;
+          int var = it->first.second;
+          expr_t d1 = it->second;
+
+          jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
+          jacobian_output << "=";
+          d1->writeOutput(jacobian_output, output_type,
+                          temp_term_union, temp_term_union_idxs, tef_terms);
+          jacobian_output << ";" << endl;
+        }
     }
 
   // 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, temp_term_empty, hessian_output, output_type, tef_terms);
-    else
-      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++)
     {
-      int eq = it->first.first;
-      int var1 = it->first.second.first;
-      int var2 = it->first.second.second;
-      expr_t d2 = it->second;
+      writeTemporaryTerms(temporary_terms_g2, temporary_terms_g2_idxs,
+                          temp_term_union, temp_term_union_idxs,
+                          hessian_tt_output, output_type, tef_terms);
+      temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
+      temp_term_union_idxs.insert(temp_term_union_idxs.end(),
+                                  temporary_terms_g2_idxs.begin(), temporary_terms_g2_idxs.end());
 
-      int id1 = getDynJacobianCol(var1);
-      int id2 = getDynJacobianCol(var2);
+      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++)
+        {
+          int eq = it->first.first;
+          int var1 = it->first.second.first;
+          int var2 = it->first.second.second;
+          expr_t d2 = it->second;
 
-      int col_nb = id1 * dynJacobianColsNbr + id2;
-      int col_nb_sym = id2 * dynJacobianColsNbr + id1;
+          int id1 = getDynJacobianCol(var1);
+          int id2 = getDynJacobianCol(var2);
 
-      ostringstream for_sym;
-      if (output_type == oJuliaDynamicModel)
-        {
-          for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
-          hessian_output << "  @inbounds " << for_sym.str() << " = ";
-          d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
-          hessian_output << endl;
-        }
-      else
-        {
-          sparseHelper(2, hessian_output, k, 0, output_type);
-          hessian_output << "=" << eq + 1 << ";" << endl;
+          int col_nb = id1 * dynJacobianColsNbr + id2;
+          int col_nb_sym = id2 * dynJacobianColsNbr + id1;
 
-          sparseHelper(2, hessian_output, k, 1, output_type);
-          hessian_output << "=" << col_nb + 1 << ";" << endl;
+          ostringstream for_sym;
+          if (output_type == oJuliaDynamicModel)
+            {
+              for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
+              hessian_output << "    @inbounds " << for_sym.str() << " = ";
+              d2->writeOutput(hessian_output, output_type, temp_term_union, temp_term_union_idxs, tef_terms);
+              hessian_output << endl;
+            }
+          else
+            {
+              sparseHelper(2, hessian_output, k, 0, output_type);
+              hessian_output << "=" << eq + 1 << ";" << endl;
 
-          sparseHelper(2, hessian_output, k, 2, output_type);
-          hessian_output << "=";
-          d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
-          hessian_output << ";" << endl;
+              sparseHelper(2, hessian_output, k, 1, output_type);
+              hessian_output << "=" << col_nb + 1 << ";" << endl;
 
-          k++;
-        }
+              sparseHelper(2, hessian_output, k, 2, output_type);
+              hessian_output << "=";
+              d2->writeOutput(hessian_output, output_type, temp_term_union, temp_term_union_idxs, tef_terms);
+              hessian_output << ";" << endl;
 
-      // Treating symetric elements
-      if (id1 != id2)
-        if (output_type == oJuliaDynamicModel)
-          hessian_output << "  @inbounds g2[" << eq + 1 << "," << col_nb_sym + 1 << "] = "
-                         << for_sym.str() << endl;
-        else
-          {
-            sparseHelper(2, hessian_output, k, 0, output_type);
-            hessian_output << "=" << eq + 1 << ";" << endl;
+              k++;
+            }
+
+          // Treating symetric elements
+          if (id1 != id2)
+            if (output_type == oJuliaDynamicModel)
+              hessian_output << "    @inbounds g2[" << eq + 1 << "," << col_nb_sym + 1 << "] = "
+                             << for_sym.str() << endl;
+            else
+              {
+                sparseHelper(2, hessian_output, k, 0, output_type);
+                hessian_output << "=" << eq + 1 << ";" << endl;
 
-            sparseHelper(2, hessian_output, k, 1, output_type);
-            hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
+                sparseHelper(2, hessian_output, k, 1, output_type);
+                hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
 
-            sparseHelper(2, hessian_output, k, 2, output_type);
-            hessian_output << "=";
-            sparseHelper(2, hessian_output, k-1, 2, output_type);
-            hessian_output << ";" << endl;
+                sparseHelper(2, hessian_output, k, 2, output_type);
+                hessian_output << "=";
+                sparseHelper(2, hessian_output, k-1, 2, output_type);
+                hessian_output << ";" << endl;
 
-            k++;
-          }
+                k++;
+              }
+        }
     }
 
   // 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, temp_term_empty, third_derivatives_output, output_type, tef_terms);
-    else
-      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++)
     {
-      int eq = it->first.first;
-      int var1 = it->first.second.first;
-      int var2 = it->first.second.second.first;
-      int var3 = it->first.second.second.second;
-      expr_t d3 = it->second;
-
-      int id1 = getDynJacobianCol(var1);
-      int id2 = getDynJacobianCol(var2);
-      int id3 = getDynJacobianCol(var3);
+      writeTemporaryTerms(temporary_terms_g3, temporary_terms_g3_idxs,
+                          temp_term_union, temp_term_union_idxs,
+                          third_derivatives_tt_output, output_type, tef_terms);
+      temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
+      temp_term_union_idxs.insert(temp_term_union_idxs.end(),
+                                  temporary_terms_g3_idxs.begin(), temporary_terms_g3_idxs.end());
 
-      // Reference column number for the g3 matrix
-      int ref_col = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
-
-      ostringstream for_sym;
-      if (output_type == oJuliaDynamicModel)
-        {
-          for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
-          third_derivatives_output << "  @inbounds " << for_sym.str() << " = ";
-          d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
-          third_derivatives_output << endl;
-        }
-      else
+      int 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++)
         {
-          sparseHelper(3, third_derivatives_output, k, 0, output_type);
-          third_derivatives_output << "=" << eq + 1 << ";" << endl;
-
-          sparseHelper(3, third_derivatives_output, k, 1, output_type);
-          third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
+          int eq = it->first.first;
+          int var1 = it->first.second.first;
+          int var2 = it->first.second.second.first;
+          int var3 = it->first.second.second.second;
+          expr_t d3 = it->second;
 
-          sparseHelper(3, third_derivatives_output, k, 2, output_type);
-          third_derivatives_output << "=";
-          d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
-          third_derivatives_output << ";" << endl;
-        }
+          int id1 = getDynJacobianCol(var1);
+          int id2 = getDynJacobianCol(var2);
+          int id3 = getDynJacobianCol(var3);
 
-      // Compute the column numbers for the 5 other permutations of (id1,id2,id3)
-      // and store them in a set (to avoid duplicates if two indexes are equal)
-      set<int> cols;
-      cols.insert(id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2);
-      cols.insert(id2 * hessianColsNbr + id1 * dynJacobianColsNbr + id3);
-      cols.insert(id2 * hessianColsNbr + id3 * dynJacobianColsNbr + id1);
-      cols.insert(id3 * hessianColsNbr + id1 * dynJacobianColsNbr + id2);
-      cols.insert(id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1);
+          // Reference column number for the g3 matrix
+          int ref_col = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
 
-      int k2 = 1; // Keeps the offset of the permutation relative to k
-      for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
-        if (*it2 != ref_col)
+          ostringstream for_sym;
           if (output_type == oJuliaDynamicModel)
-            third_derivatives_output << "  @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
-                                     << for_sym.str() << endl;
+            {
+              for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
+              third_derivatives_output << "    @inbounds " << for_sym.str() << " = ";
+              d3->writeOutput(third_derivatives_output, output_type, temp_term_union, temp_term_union_idxs, tef_terms);
+              third_derivatives_output << endl;
+            }
           else
             {
-              sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
+              sparseHelper(3, third_derivatives_output, k, 0, output_type);
               third_derivatives_output << "=" << eq + 1 << ";" << endl;
 
-              sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
-              third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
+              sparseHelper(3, third_derivatives_output, k, 1, output_type);
+              third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
 
-              sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
-              third_derivatives_output << "=";
               sparseHelper(3, third_derivatives_output, k, 2, output_type);
+              third_derivatives_output << "=";
+              d3->writeOutput(third_derivatives_output, output_type, temp_term_union, temp_term_union_idxs, tef_terms);
               third_derivatives_output << ";" << endl;
-
-              k2++;
             }
-      k += k2;
+
+          // Compute the column numbers for the 5 other permutations of (id1,id2,id3)
+          // and store them in a set (to avoid duplicates if two indexes are equal)
+          set<int> cols;
+          cols.insert(id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2);
+          cols.insert(id2 * hessianColsNbr + id1 * dynJacobianColsNbr + id3);
+          cols.insert(id2 * hessianColsNbr + id3 * dynJacobianColsNbr + id1);
+          cols.insert(id3 * hessianColsNbr + id1 * dynJacobianColsNbr + id2);
+          cols.insert(id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1);
+
+          int k2 = 1; // Keeps the offset of the permutation relative to k
+          for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
+            if (*it2 != ref_col)
+              if (output_type == oJuliaDynamicModel)
+                third_derivatives_output << "    @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
+                                         << for_sym.str() << endl;
+              else
+                {
+                  sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
+                  third_derivatives_output << "=" << eq + 1 << ";" << endl;
+
+                  sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
+                  third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
+
+                  sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
+                  third_derivatives_output << "=";
+                  sparseHelper(3, third_derivatives_output, k, 2, output_type);
+                  third_derivatives_output << ";" << endl;
+
+                  k2++;
+                }
+          k += k2;
+        }
     }
 
   if (output_type == oMatlabDynamicModel)
@@ -2403,192 +2517,306 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
       map<string, string> tmp_paren_vars;
       bool message_printed = false;
       fixNestedParenthesis(model_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(model_local_vars_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(model_tt_output, tmp_paren_vars, message_printed);
       fixNestedParenthesis(jacobian_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(jacobian_tt_output, tmp_paren_vars, message_printed);
       fixNestedParenthesis(hessian_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(hessian_tt_output, tmp_paren_vars, message_printed);
       fixNestedParenthesis(third_derivatives_output, tmp_paren_vars, message_printed);
-
-      DynamicOutput << "%" << endl
-                    << "% Model equations" << endl
-                    << "%" << endl
-                    << endl;
-      writeVarExpectationCalls(DynamicOutput);
-      DynamicOutput << "residual = zeros(" << nrows << ", 1);" << endl
-                    << model_local_vars_output.str()
-                    << model_output.str()
-        // Writing initialization instruction for matrix g1
-                    << "if nargout >= 2," << endl
-                    << "  g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");" << endl
-                    << endl
-                    << "  %" << endl
-                    << "  % Jacobian matrix" << endl
-                    << "  %" << endl
-                    << endl
-                    << jacobian_output.str()
-                    << endl
-
-        // Initialize g2 matrix
-                    << "if nargout >= 3," << endl
-                    << "  %" << endl
-                    << "  % Hessian matrix" << endl
-                    << "  %" << endl
-                    << endl;
+      fixNestedParenthesis(third_derivatives_tt_output, tmp_paren_vars, message_printed);
+
+      ostringstream init_output, end_output;
+      init_output << "residual = zeros(" << nrows << ", 1);";
+      writeDynamicModelHelper(dynamic_basename + "_resid", "residual",
+                              dynamic_basename + "_resid_tt",
+                              temporary_terms_res_idxs.size(),
+                              "", init_output, end_output,
+                              model_output, model_tt_output);
+
+      init_output.str(string());
+      init_output.clear();
+      init_output << "g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");";
+      writeDynamicModelHelper(dynamic_basename + "_g1", "g1",
+                              dynamic_basename + "_g1_tt",
+                              temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size(),
+                              dynamic_basename + "_resid_tt",
+                              init_output, end_output,
+                              jacobian_output, jacobian_tt_output);
+      writeWrapperFunctions(dynamic_basename, "g1");
+
+      init_output.str(string());
+      init_output.clear();
       if (second_derivatives.size())
-        DynamicOutput << "  v2 = zeros(" << NNZDerivatives[1] << ",3);" << endl
-                      << hessian_output.str()
-                      << "  g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << nrows << "," << hessianColsNbr << ");" << endl;
-      else // Either hessian is all zero, or we didn't compute it
-        DynamicOutput << "  g2 = sparse([],[],[]," << nrows << "," << hessianColsNbr << ");" << endl;
-
-      // Initialize g3 matrix
-      DynamicOutput << "if nargout >= 4," << endl
-                    << "  %" << endl
-                    << "  % Third order derivatives" << endl
-                    << "  %" << endl
-                    << endl;
+        {
+          init_output << "v2 = zeros(" << NNZDerivatives[1] << ",3);";
+          end_output << "g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << nrows << "," << hessianColsNbr << ");";
+        }
+      else
+        init_output << "g2 = sparse([],[],[]," << nrows << "," << hessianColsNbr << ");";
+      writeDynamicModelHelper(dynamic_basename + "_g2", "g2",
+                              dynamic_basename + "_g2_tt",
+                              temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size()
+                              + temporary_terms_g2_idxs.size(),
+                              dynamic_basename + "_g1_tt",
+                              init_output, end_output,
+                              hessian_output, hessian_tt_output);
+      writeWrapperFunctions(dynamic_basename, "g2");
+
+      init_output.str(string());
+      init_output.clear();
+      end_output.str(string());
+      end_output.clear();
       int ncols = hessianColsNbr * dynJacobianColsNbr;
       if (third_derivatives.size())
-        DynamicOutput << "  v3 = zeros(" << NNZDerivatives[2] << ",3);" << endl
-                      << third_derivatives_output.str()
-                      << "  g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");" << endl;
-      else // Either 3rd derivatives is all zero, or we didn't compute it
-        DynamicOutput << "  g3 = sparse([],[],[]," << nrows << "," << ncols << ");" << endl;
+        {
+          init_output << "v3 = zeros(" << NNZDerivatives[2] << ",3);";
+          end_output << "g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");";
+        }
+      else
+        init_output << "g3 = sparse([],[],[]," << nrows << "," << ncols << ");";
+      writeDynamicModelHelper(dynamic_basename + "_g3", "g3",
+                              dynamic_basename + "_g3_tt",
+                              temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size()
+                              + temporary_terms_g2_idxs.size() + temporary_terms_g3_idxs.size(),
+                              dynamic_basename + "_g2_tt",
+                              init_output, end_output,
+                              third_derivatives_output, third_derivatives_tt_output);
+      writeWrapperFunctions(dynamic_basename, "g3");
+    }
+        /*
 
-      DynamicOutput << "end" << endl
-                    << "end" << endl
-                    << "end" << endl;
-    }
   else if (output_type == oCDynamicModel)
     {
       DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, double *g1, double *v2, double *v3)" << endl
                     << "{" << endl
                     << "  double lhs, rhs;" << endl
                     << endl
-                    << "  /* Residual equations */" << endl
+                    << "  * Residual equations *" << endl
                     << model_local_vars_output.str()
+                    << model_tt_output.str()
                     << model_output.str()
-                    << "  /* Jacobian  */" << endl
+                    << "  * Jacobian  *" << endl
                     << "  if (g1 == NULL)" << endl
                     << "    return;" << endl
                     << endl
+                    << jacobian_tt_output.str()
                     << jacobian_output.str()
                     << endl;
 
       if (second_derivatives.size())
-        DynamicOutput << "  /* Hessian for endogenous and exogenous variables */" << endl
+        DynamicOutput << "  * Hessian for endogenous and exogenous variables *" << endl
                       << "  if (v2 == NULL)" << endl
                       << "    return;" << endl
                       << endl
+                      << hessian_tt_output.str()
                       << hessian_output.str()
                       << endl;
 
       if (third_derivatives.size())
-        DynamicOutput << "  /* Third derivatives for endogenous and exogenous variables */" << endl
+        DynamicOutput << "  * Third derivatives for endogenous and exogenous variables *" << endl
                       << "  if (v3 == NULL)" << endl
                       << "    return;" << endl
                       << endl
+                      << third_derivatives_tt_output.str()
                       << third_derivatives_output.str()
                       << endl;
 
       DynamicOutput << "}" << endl << endl;
     }
+        */
   else
     {
-      ostringstream comments;
-      comments << "## Function Arguments" << endl
-               << endl
-               << "## Input" << endl
-               << " 1 y:            Array{Float64, num_dynamic_vars, 1}             Vector of endogenous variables in the order stored" << endl
-               << "                                                                 in model_.lead_lag_incidence; see the manual" << endl
-               << " 2 x:            Array{Float64, nperiods, length(model_.exo)}    Matrix of exogenous variables (in declaration order)" << endl
-               << "                                                                 for all simulation periods" << endl
-               << " 3 params:       Array{Float64, length(model_.param), 1}         Vector of parameter values in declaration order" << endl
-               << " 4 steady_state:" << endl
-               << " 5 it_:          Int                                             Time period for exogenous variables for which to evaluate the model" << endl
-               << endl
-               << "## Output" << endl
-               << " 6 residual:     Array(Float64, model_.eq_nbr, 1)                Vector of residuals of the dynamic model equations in" << endl
-               << "                                                                 order of declaration of the equations." << endl;
-
-      DynamicOutput << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, "
-                    << "params::Vector{Float64}," << endl
-                    << "                  steady_state::Vector{Float64}, it_::Int, "
-                    << "residual::Vector{Float64})" << endl
-                    << "#=" << endl << comments.str() << "=#" << endl
-                    << "  @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
-                    << "  @assert length(params) == " << symbol_table.param_nbr() << endl
-                    << "  @assert length(residual) == " << nrows << endl
-                    << "  #" << endl
-                    << "  # Model equations" << endl
-                    << "  #" << endl
-                    << model_local_vars_output.str()
-                    << model_output.str()
-                    << "end" << endl << endl
-                    << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, "
-                    << "params::Vector{Float64}," << endl
-                    << "                  steady_state::Vector{Float64}, it_::Int, "
-                    << "residual::Vector{Float64}," << endl
-                    << "                  g1::Matrix{Float64})" << endl;
-
-      comments << " 7 g1:           Array(Float64, model_.eq_nbr, num_dynamic_vars) Jacobian matrix of the dynamic model equations;" << endl
-               << "                                                                 rows: equations in order of declaration" << endl
-               << "                                                                 columns: variables in order stored in model_.lead_lag_incidence" << endl;
-
-      DynamicOutput << "#=" << endl << comments.str() << "=#" << endl
-                    << "  @assert size(g1) == (" << nrows << ", " << dynJacobianColsNbr << ")" << endl
-                    << "  fill!(g1, 0.0)" << endl
-                    << "  dynamic!(y, x, params, steady_state, it_, residual)" << endl
-                    << model_local_vars_output.str()
-                    << "  #" << endl
-                    << "  # Jacobian matrix" << endl
-                    << "  #" << endl
-                    << jacobian_output.str()
-                    << "end" << endl << endl
-                    << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, "
-                    << "params::Vector{Float64}," << endl
-                    << "                  steady_state::Vector{Float64}, it_::Int, "
-                    << "residual::Vector{Float64}," << endl
-                    << "                  g1::Matrix{Float64}, g2::Matrix{Float64})" << endl;
-
-      comments << " 8 g2:           spzeros(model_.eq_nbr, (num_dynamic_vars)^2)    Hessian matrix of the dynamic model equations;" << endl
-               << "                                                                 rows: equations in order of declaration" << endl
-               << "                                                                 columns: variables in order stored in model_.lead_lag_incidence" << endl;
-
-      DynamicOutput << "#=" << endl << comments.str() << "=#" << endl
-                    << "  @assert size(g2) == (" << nrows << ", " << hessianColsNbr << ")" << endl
-                    << "  fill!(g2, 0.0)" << endl
-                    << "  dynamic!(y, x, params, steady_state, it_, residual, g1)" << endl;
-      if (second_derivatives.size())
-        DynamicOutput << model_local_vars_output.str()
-                      << "  #" << endl
-                      << "  # Hessian matrix" << endl
-                      << "  #" << endl
-                      << hessian_output.str();
+      string filename =  dynamic_basename + "Dynamic.jl";
+      ofstream output;
+      output.open(filename.c_str(), ios::out | ios::binary);
+      if (!output.is_open())
+        {
+          cerr << "Error: Can't open file " << filename << " for writing" << endl;
+          exit(EXIT_FAILURE);
+        }
 
-      // Initialize g3 matrix
+      output << "module " << dynamic_basename << "Dynamic" << endl
+             << "#" << endl
+             << "# NB: this file was automatically generated by Dynare" << endl
+             << "#     from " << dynamic_basename << ".mod" << endl
+             << "#" << endl
+             << "using Utils" << endl << endl
+             << "export dynamic!, dynamicResid!, dynamicG1!, dynamicG2!, dynamicG3!" << endl << endl
+             << "#=" << endl
+             << "# The comments below apply to all functions contained in this module #" << endl
+             << "  NB: The arguments contained on the first line of the function" << endl
+             << "      definition are those that are modified in place" << endl << endl
+             << "## Exported Functions ##" << endl
+             << "  dynamic!      : Wrapper function; computes residuals, Jacobian, Hessian," << endl
+             << "                  and third derivatives depending on the arguments provided" << endl
+             << "  dynamicResid! : Computes the dynamic model residuals" << endl
+             << "  dynamicG1!    : Computes the dynamic model Jacobian" << endl
+             << "  dynamicG2!    : Computes the dynamic model Hessian" << endl
+             << "  dynamicG3!    : Computes the dynamic model third derivatives" << endl << endl
+             << "## Local Functions ##" << endl
+             << "  dynamicResidTT! : Computes the dynamic model temporary terms for the residuals" << endl
+             << "  dynamicG1TT!    : Computes the dynamic model temporary terms for the Jacobian" << endl
+             << "  dynamicG2TT!    : Computes the dynamic model temporary terms for the Hessian" << endl
+             << "  dynamicG3TT!    : Computes the dynamic model temporary terms for the third derivatives" << endl << endl
+             << "## Function Arguments ##" << endl
+             << "  T            : Vector{Float64, num_temp_terms, 1}                Vector of Temporary Terms" << endl
+             << "  y            : Vector{Float64, num_dynamic_vars, 1}              Vector of endogenous variables in the order stored" << endl
+             << "                                                                   in model_.lead_lag_incidence; see the manual" << endl
+             << "  x            : Vector{Float64, nperiods, length(model_.exo)}     Matrix of exogenous variables (in declaration order)" << endl
+             << "                                                                   for all simulation periods" << endl
+             << "  params       : Vector{Float64, length(model_.param), 1}          Vector of parameter values in declaration order" << endl
+             << "  steady_state :" << endl
+             << "  it_          : Int                                               Time period for exogenous variables for which to evaluate the model" << endl
+             << "  residual     : Vector{Float64, model_.eq_nbr, 1}                 Vector of residuals of the dynamic model equations in" << endl
+             << "                                                                   order of declaration of the equations." << endl
+             << "  g1           : Vector{Float64, model_.eq_nbr, num_dynamic_vars}  Jacobian matrix of the dynamic model equations" << endl
+             << "                                                                   rows: equations in order of declaration" << endl
+             << "                                                                   columns: variables in order stored in model_.lead_lag_incidence" << endl
+             << "  g2           : spzeros(model_.eq_nbr, (num_dynamic_vars)^2)      Hessian matrix of the dynamic model equations" << endl
+             << "                                                                   rows: equations in order of declaration" << endl
+             << "                                                                   columns: variables in order stored in model_.lead_lag_incidence" << endl
+             << "  g3           : spzeros(model_.eq_nbr, (num_dynamic_vars)^3)      Third order derivative matrix of the dynamic model equations;" << endl
+             << "                                                                   rows: equations in order of declaration" << endl
+             << "                                                                   columns: variables in order stored in model_.lead_lag_incidence" << endl
+             << "=#" << endl << endl;
+
+      // dynamicResidTT!
+      output << "function dynamicResidTT!(T::Vector{Float64}," << endl
+             << "                         y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+             << model_tt_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // dynamic!
+      output << "function dynamicResid!(T::Vector{Float64}, residual::Vector{Float64}," << endl
+             << "                       y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
+             << "    @assert length(T) >= " << temporary_terms_res_idxs.size() << endl
+             << "    @assert length(residual) == " << nrows << endl
+             << "    @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
+             << "    @assert length(params) == " << symbol_table.param_nbr() << endl
+             << "    if T_flag" << endl
+             << "        dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl
+             << "    end" << endl
+             << model_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // dynamicG1TT!
+      output << "function dynamicG1TT!(T::Vector{Float64}," << endl
+             << "                      y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+             << "    dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl
+             << jacobian_tt_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // dynamicG1!
+      output << "function dynamicG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl
+             << "                    y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
+             << "    @assert length(T) >= "
+             << temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size() << endl
+             << "    @assert size(g1) == (" << nrows << ", " << dynJacobianColsNbr << ")" << endl
+             << "    @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
+             << "    @assert length(params) == " << symbol_table.param_nbr() << endl
+             << "    if T_flag" << endl
+             << "        dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl
+             << "    end" << endl
+             << "    fill!(g1, 0.0)" << endl
+             << jacobian_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // dynamicG2TT!
+      output << "function dynamicG2TT!(T::Vector{Float64}," << endl
+             << "                      y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+             << "    dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl
+             << hessian_tt_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // dynamicG2!
+      output << "function dynamicG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl
+             << "                    y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
+             << "    @assert length(T) >= " << temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size() + temporary_terms_g2_idxs.size() << endl
+             << "    @assert size(g2) == (" << nrows << ", " << hessianColsNbr << ")" << endl
+             << "    @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
+             << "    @assert length(params) == " << symbol_table.param_nbr() << endl
+             << "    if T_flag" << endl
+             << "        dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl
+             << "    end" << endl
+             << "    fill!(g2, 0.0)" << endl
+             << hessian_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // dynamicG3TT!
+      output << "function dynamicG3TT!(T::Vector{Float64}," << endl
+             << "                      y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+             << "    dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl
+             << third_derivatives_tt_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // dynamicG3!
       int ncols = hessianColsNbr * dynJacobianColsNbr;
-      DynamicOutput << "end" << endl << endl
-                    << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, "
-                    << "params::Vector{Float64}," << endl
-                    << "                  steady_state::Vector{Float64}, it_::Int, "
-                    << "residual::Vector{Float64}," << endl
-                    << "                  g1::Matrix{Float64}, g2::Matrix{Float64}, g3::Matrix{Float64})" << endl;
-
-      comments << " 9 g3:           spzeros(model_.eq_nbr, (num_dynamic_vars)^3)    Third order derivative matrix of the dynamic model equations;" << endl
-               << "                                                                 rows: equations in order of declaration" << endl
-               << "                                                                 columns: variables in order stored in model_.lead_lag_incidence" << endl;
-
-      DynamicOutput << "#=" << endl << comments.str() << "=#" << endl
-                    << "  @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
-                    << "  fill!(g3, 0.0)" << endl
-                    << "  dynamic!(y, x, params, steady_state, it_, residual, g1, g2)" << endl;
-      if (third_derivatives.size())
-        DynamicOutput << model_local_vars_output.str()
-                      << "  #" << endl
-                      << "  # Third order derivatives" << endl
-                      << "  #" << endl
-                      << third_derivatives_output.str();
-      DynamicOutput << "end" << endl;
+      output << "function dynamicG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl
+             << "                    y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
+             << "    @assert length(T) >= "
+             << temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size() + temporary_terms_g2_idxs.size() + temporary_terms_g3_idxs.size() << endl
+             << "    @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
+             << "    @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl
+             << "    @assert length(params) == " << symbol_table.param_nbr() << endl
+             << "    if T_flag" << endl
+             << "      dynamicG3TT!(T, y, x, params, steady_state, it_)" << endl
+             << "    end" << endl
+             << "    fill!(g3, 0.0)" << endl
+             << third_derivatives_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // dynamic!
+      output << "function dynamic!(T::Vector{Float64}, residual::Vector{Float64}," << endl
+             << "                  y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+             << "    dynamicResid!(T, residual, y, x, params, steady_state, it_, true)" << endl
+             << "    return nothing" << endl
+             << "end" << endl
+             << endl
+             << "function dynamic!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}," << endl
+             << "                  y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+             << "    dynamicG1!(T, g1, y, x, params, steady_state, it_, true)" << endl
+             << "    dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl
+             << "    return nothing" << endl
+             << "end" << endl
+             << endl
+             << "function dynamic!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl
+             << "                  y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+             << "    dynamicG2!(T, g2, y, x, params, steady_state, it_, true)" << endl
+             << "    dynamicG1!(T, g1, y, x, params, steady_state, it_, false)" << endl
+             << "    dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl
+             << "    return nothing" << endl
+             << "end" << endl
+             << endl
+             << "function dynamic!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}, g3::Matrix{Float64}," << endl
+             << "                  y::Vector{Float64}, x::Matrix{Float64}, "
+             << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+             << "    dynamicG3!(T, g3, y, x, params, steady_state, it_, true)" << endl
+             << "    dynamicG2!(T, g2, y, x, params, steady_state, it_, false)" << endl
+             << "    dynamicG1!(T, g1, y, x, params, steady_state, it_, false)" << endl
+             << "    dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl
+             << "    return nothing" << endl
+             << "end" << endl
+             << "end" << endl;
+      output.close();
     }
 }
 
@@ -4031,13 +4259,16 @@ DynamicModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputTyp
 {
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temporary_terms;
+  temporary_terms_idxs_t temporary_terms_idxs;
   for (int i = 0; i < (int) aux_equations.size(); i++)
     if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
       dynamic_cast<ExprNode *>(aux_equations[i])->writeExternalFunctionOutput(output, output_type,
-                                                                              temporary_terms, tef_terms);
+                                                                              temporary_terms,
+                                                                              temporary_terms_idxs,
+                                                                              tef_terms);
   for (int i = 0; i < (int) aux_equations.size(); i++)
     {
-      dynamic_cast<ExprNode *>(aux_equations[i])->writeOutput(output, output_type, temporary_terms, tef_terms);
+      dynamic_cast<ExprNode *>(aux_equations[i])->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ";" << endl;
     }
 }
@@ -5266,11 +5497,11 @@ DynamicModel::isChecksumMatching(const string &basename) const
       if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
         {
           buffer << "lhs =";
-          lhs->writeOutput(buffer, buffer_type, temporary_terms);
+          lhs->writeOutput(buffer, buffer_type, temporary_terms, temporary_terms_idxs);
           buffer << ";" << endl;
 
           buffer << "rhs =";
-          rhs->writeOutput(buffer, buffer_type, temporary_terms);
+          rhs->writeOutput(buffer, buffer_type, temporary_terms, temporary_terms_idxs);
           buffer << ";" << endl;
 
           buffer << "residual" << LEFT_ARRAY_SUBSCRIPT(buffer_type)
@@ -5284,7 +5515,7 @@ DynamicModel::isChecksumMatching(const string &basename) const
                  << eq + ARRAY_SUBSCRIPT_OFFSET(buffer_type)
                  << RIGHT_ARRAY_SUBSCRIPT(buffer_type)
                  << " = ";
-          lhs->writeOutput(buffer, buffer_type, temporary_terms);
+          lhs->writeOutput(buffer, buffer_type, temporary_terms, temporary_terms_idxs);
           buffer << ";" << endl;
         }
     }
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 7ff8475f4b83b0d86fbc0d054c5cc639ba1f3e9b..591eae8d5d8699d7cd2387c90d5a817acd9c85db 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -109,6 +109,8 @@ private:
   //! Writes the dynamic model equations and its derivatives
   /*! \todo add third derivatives handling in C output */
   void writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const;
+  void writeDynamicModel(const string &dynamic_basename, bool use_dll, bool julia) const;
+  void writeDynamicModel(const string &dynamic_basename, ostream &DynamicOutput, bool use_dll, bool julia) const;
   //! Writes the Block reordred structure of the model in M output
   void writeModelEquationsOrdered_M(const string &dynamic_basename) const;
   //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
@@ -235,6 +237,14 @@ private:
   //!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
   vector<pair<int, int> > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
 
+  void writeWrapperFunctions(const string &name, const string &ending) const;
+  void writeDynamicModelHelper(const string &name, const string &retvalname,
+                               const string &name_tt, size_t ttlen,
+                               const string &previous_tt_name,
+                               const ostringstream &init_s,
+                               const ostringstream &end_s,
+                               const ostringstream &s, const ostringstream &s_tt) const;
+
 public:
   DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_argx);
   //! Adds a variable node
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index f2340b6508b48a5c7e6f95bb6a8611f68b8f2155..4b84103bbcfa08cdb45a01e4523429cc0fbefef6 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -171,8 +171,16 @@ ExprNode::writeOutput(ostream &output, ExprNodeOutputType output_type) const
 void
 ExprNode::writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const
 {
+  temporary_terms_idxs_t temporary_terms_idxs;
   deriv_node_temp_terms_t tef_terms;
-  writeOutput(output, output_type, temporary_terms, tef_terms);
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
+void
+ExprNode::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;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 }
 
 void
@@ -187,6 +195,7 @@ ExprNode::compile(ostream &CompileCode, unsigned int &instruction_number,
 void
 ExprNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                       const temporary_terms_t &temporary_terms,
+                                      const temporary_terms_idxs_t &temporary_terms_idxs,
                                       deriv_node_temp_terms_t &tef_terms) const
 {
   // Nothing to do
@@ -344,13 +353,28 @@ void
 NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                           const temporary_terms_t &temporary_terms,
                           deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
+void
+NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                          const temporary_terms_t &temporary_terms,
+                          const temporary_terms_idxs_t &temporary_terms_idxs,
+                          deriv_node_temp_terms_t &tef_terms) const
 {
   temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<NumConstNode *>(this));
   if (it != temporary_terms.end())
     if (output_type == oMatlabDynamicModelSparse)
       output << "T" << idx << "(it_)";
     else
-      output << "T" << idx;
+      if (temporary_terms_idxs.empty() || IS_C(output_type))
+        output << "T" << idx;
+      else
+        output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+               << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+               << RIGHT_ARRAY_SUBSCRIPT(output_type);
   else
     output << datatree.num_constants.get(id);
 }
@@ -391,6 +415,11 @@ NumConstNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int> >
 {
 }
 
+void
+NumConstNode::collectModelLocalVariables(map<int, expr_t> &result) const
+{
+}
+
 pair<int, expr_t >
 NumConstNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
@@ -739,6 +768,7 @@ VariableNode::writeJsonOutput(ostream &output,
 void
 VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                           const temporary_terms_t &temporary_terms,
+                          const temporary_terms_idxs_t &temporary_terms_idxs,
                           deriv_node_temp_terms_t &tef_terms) const
 {
   // If node is a temporary term
@@ -748,7 +778,14 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       if (output_type == oMatlabDynamicModelSparse)
         output << "T" << idx << "(it_)";
       else
-        output << "T" << idx;
+        {
+        if (temporary_terms_idxs.empty() || IS_C(output_type))
+          output << "T" << idx;
+        else
+          output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                 << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type);
+        }
       return;
     }
 
@@ -790,7 +827,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           || output_type == oCDynamicSteadyStateOperator)
         {
           output << "(";
-          datatree.local_variables_table[symb_id]->writeOutput(output, output_type, temporary_terms, tef_terms);
+          datatree.local_variables_table[symb_id]->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ")";
         }
       else
@@ -986,6 +1023,15 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     }
 }
 
+void
+VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                          const temporary_terms_t &temporary_terms,
+                          deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 expr_t
 VariableNode::substituteStaticAuxiliaryVariable() const
 {
@@ -1107,6 +1153,16 @@ VariableNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int> >
     datatree.local_variables_table[symb_id]->collectDynamicVariables(type_arg, result);
 }
 
+void
+VariableNode::collectModelLocalVariables(map<int, expr_t> &result) const
+{
+  if (type == eModelLocalVariable)
+    {
+      result[symb_id] = const_cast<VariableNode *>(this);
+      datatree.local_variables_table[symb_id]->collectModelLocalVariables(result);
+    }
+}
+
 pair<int, expr_t>
 VariableNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
@@ -2249,6 +2305,7 @@ UnaryOpNode::writeJsonOutput(ostream &output,
 void
 UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                          const temporary_terms_t &temporary_terms,
+                         const temporary_terms_idxs_t &temporary_terms_idxs,
                          deriv_node_temp_terms_t &tef_terms) const
 {
   // If node is a temporary term
@@ -2258,7 +2315,12 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       if (output_type == oMatlabDynamicModelSparse)
         output << "T" << idx << "(it_)";
       else
-        output << "T" << idx;
+        if (temporary_terms_idxs.empty() || IS_C(output_type))
+          output << "T" << idx;
+        else
+          output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                 << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type);
       return;
     }
 
@@ -2358,7 +2420,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           break;
         }
       output << "(";
-      arg->writeOutput(output, new_output_type, temporary_terms, tef_terms);
+      arg->writeOutput(output, new_output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ")";
       return;
     case oSteadyStateParamDeriv:
@@ -2431,7 +2493,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     }
 
   // Write argument
-  arg->writeOutput(output, output_type, temporary_terms, tef_terms);
+  arg->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 
   if (close_parenthesis)
     output << RIGHT_PAR(output_type);
@@ -2441,12 +2503,22 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     output << RIGHT_PAR(output_type);
 }
 
+void
+UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                         const temporary_terms_t &temporary_terms,
+                         deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 void
 UnaryOpNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                          const temporary_terms_t &temporary_terms,
+                                         const temporary_terms_idxs_t &temporary_terms_idxs,
                                          deriv_node_temp_terms_t &tef_terms) const
 {
-  arg->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
+  arg->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 }
 
 void
@@ -2579,6 +2651,12 @@ UnaryOpNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &
     arg->collectDynamicVariables(type_arg, result);
 }
 
+void
+UnaryOpNode::collectModelLocalVariables(map<int, expr_t> &result) const
+{
+  arg->collectModelLocalVariables(result);
+}
+
 pair<int, expr_t>
 UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
@@ -3857,6 +3935,7 @@ BinaryOpNode::writeJsonOutput(ostream &output,
 void
 BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                           const temporary_terms_t &temporary_terms,
+                          const temporary_terms_idxs_t &temporary_terms_idxs,
                           deriv_node_temp_terms_t &tef_terms) const
 {
   // If current node is a temporary term
@@ -3866,7 +3945,12 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       if (output_type == oMatlabDynamicModelSparse)
         output << "T" << idx << "(it_)";
       else
-        output << "T" << idx;
+        if (temporary_terms_idxs.empty() || IS_C(output_type))
+          output << "T" << idx;
+        else
+          output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                 << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type);
       return;
     }
 
@@ -3874,16 +3958,16 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
   if (op_code == oPowerDeriv)
     {
       if (IS_LATEX(output_type))
-        unpackPowerDeriv()->writeOutput(output, output_type, temporary_terms, tef_terms);
+        unpackPowerDeriv()->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       else
         {
           if (output_type == oJuliaStaticModel || output_type == oJuliaDynamicModel)
             output << "get_power_deriv(";
           else
             output << "getPowerDeriv(";
-          arg1->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ",";
-          arg2->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << "," << powerDerivOrder << ")";
         }
       return;
@@ -3906,9 +3990,9 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         default:
           ;
         }
-      arg1->writeOutput(output, output_type, temporary_terms, tef_terms);
+      arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ",";
-      arg2->writeOutput(output, output_type, temporary_terms, tef_terms);
+      arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ")";
       return;
     }
@@ -3932,7 +4016,7 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     }
 
   // Write left argument
-  arg1->writeOutput(output, output_type, temporary_terms, tef_terms);
+  arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 
   if (close_parenthesis)
     output << RIGHT_PAR(output_type);
@@ -4025,7 +4109,7 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     }
 
   // Write right argument
-  arg2->writeOutput(output, output_type, temporary_terms, tef_terms);
+  arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 
   if (IS_LATEX(output_type) && (op_code == oPower || op_code == oDivide))
     output << "}";
@@ -4034,13 +4118,23 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     output << RIGHT_PAR(output_type);
 }
 
+void
+BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                         const temporary_terms_t &temporary_terms,
+                         deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 void
 BinaryOpNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                           const temporary_terms_t &temporary_terms,
+                                          const temporary_terms_idxs_t &temporary_terms_idxs,
                                           deriv_node_temp_terms_t &tef_terms) const
 {
-  arg1->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
-  arg2->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
+  arg1->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+  arg2->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 }
 
 void
@@ -4072,6 +4166,13 @@ BinaryOpNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int> >
   arg2->collectDynamicVariables(type_arg, result);
 }
 
+void
+BinaryOpNode::collectModelLocalVariables(map<int, expr_t> &result) const
+{
+  arg1->collectModelLocalVariables(result);
+  arg2->collectModelLocalVariables(result);
+}
+
 expr_t
 BinaryOpNode::Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const
 {
@@ -5188,13 +5289,22 @@ TrinaryOpNode::writeJsonOutput(ostream &output,
 void
 TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                            const temporary_terms_t &temporary_terms,
+                           const temporary_terms_idxs_t &temporary_terms_idxs,
                            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));
   if (it != temporary_terms.end())
     {
-      output << "T" << idx;
+      if (temporary_terms_idxs.empty())
+        output << "T" << idx;
+      else
+        if (temporary_terms_idxs.empty() || IS_C(output_type))
+          output << "T" << idx;
+        else
+          output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                 << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type);
       return;
     }
 
@@ -5205,21 +5315,21 @@ TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         {
           // In C, there is no normcdf() primitive, so use erf()
           output << "(0.5*(1+erf(((";
-          arg1->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ")-(";
-          arg2->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << "))/(";
-          arg3->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg3->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ")/M_SQRT2)))";
         }
       else
         {
           output << "normcdf(";
-          arg1->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ",";
-          arg2->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ",";
-          arg3->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg3->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ")";
         }
       break;
@@ -5228,37 +5338,47 @@ TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         {
           //(1/(v3*sqrt(2*M_PI)*exp(pow((v1-v2)/v3,2)/2)))
           output << "(1/(";
-          arg3->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg3->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << "*sqrt(2*M_PI)*exp(pow((";
-          arg1->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << "-";
-          arg2->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ")/";
-          arg3->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg3->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ",2)/2)))";
         }
       else
         {
           output << "normpdf(";
-          arg1->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ",";
-          arg2->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ",";
-          arg3->writeOutput(output, output_type, temporary_terms, tef_terms);
+          arg3->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ")";
         }
       break;
     }
 }
 
+void
+TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                           const temporary_terms_t &temporary_terms,
+                           deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 void
 TrinaryOpNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs,
                                            deriv_node_temp_terms_t &tef_terms) const
 {
-  arg1->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
-  arg2->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
-  arg3->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
+  arg1->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+  arg2->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+  arg3->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 }
 
 void
@@ -5294,6 +5414,14 @@ TrinaryOpNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int> >
   arg3->collectDynamicVariables(type_arg, result);
 }
 
+void
+TrinaryOpNode::collectModelLocalVariables(map<int, expr_t> &result) const
+{
+  arg1->collectModelLocalVariables(result);
+  arg2->collectModelLocalVariables(result);
+  arg3->collectModelLocalVariables(result);
+}
+
 pair<int, expr_t>
 TrinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
@@ -5709,6 +5837,14 @@ AbstractExternalFunctionNode::collectDynamicVariables(SymbolType type_arg, set<p
     (*it)->collectDynamicVariables(type_arg, result);
 }
 
+void
+AbstractExternalFunctionNode::collectModelLocalVariables(map<int, expr_t> &result) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin();
+       it != arguments.end(); it++)
+    (*it)->collectModelLocalVariables(result);
+}
+
 void
 AbstractExternalFunctionNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
 {
@@ -6055,6 +6191,7 @@ AbstractExternalFunctionNode::normalizeEquation(int var_endo, vector<pair<int, p
 void
 AbstractExternalFunctionNode::writeExternalFunctionArguments(ostream &output, ExprNodeOutputType output_type,
                                                              const temporary_terms_t &temporary_terms,
+                                                             const temporary_terms_idxs_t &temporary_terms_idxs,
                                                              deriv_node_temp_terms_t &tef_terms) const
 {
   for (vector<expr_t>::const_iterator it = arguments.begin();
@@ -6063,7 +6200,7 @@ AbstractExternalFunctionNode::writeExternalFunctionArguments(ostream &output, Ex
       if (it != arguments.begin())
         output << ",";
 
-      (*it)->writeOutput(output, output_type, temporary_terms, tef_terms);
+      (*it)->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
     }
 }
 
@@ -6267,6 +6404,7 @@ ExternalFunctionNode::writeJsonOutput(ostream &output,
 void
 ExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                                   const temporary_terms_t &temporary_terms,
+                                  const temporary_terms_idxs_t &temporary_terms_idxs,
                                   deriv_node_temp_terms_t &tef_terms) const
 {
   if (output_type == oMatlabOutsideModel || output_type == oSteadyStateFile
@@ -6276,7 +6414,7 @@ ExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_typ
       string name = IS_LATEX(output_type) ? datatree.symbol_table.getTeXName(symb_id)
         : datatree.symbol_table.getName(symb_id);
       output << name << "(";
-      writeExternalFunctionArguments(output, output_type, temporary_terms, tef_terms);
+      writeExternalFunctionArguments(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ")";
       return;
     }
@@ -6287,7 +6425,12 @@ ExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_typ
       if (output_type == oMatlabDynamicModelSparse)
         output << "T" << idx << "(it_)";
       else
-        output << "T" << idx;
+        if (temporary_terms_idxs.empty() || IS_C(output_type))
+          output << "T" << idx;
+        else
+          output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                 << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type);
       return;
     }
 
@@ -6296,9 +6439,19 @@ ExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_typ
   output << "TEF_" << getIndxInTefTerms(symb_id, tef_terms);
 }
 
+void
+ExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                         const temporary_terms_t &temporary_terms,
+                         deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 void
 ExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                                   const temporary_terms_t &temporary_terms,
+                                                  const temporary_terms_idxs_t &temporary_terms_idxs,
                                                   deriv_node_temp_terms_t &tef_terms) const
 {
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
@@ -6306,7 +6459,7 @@ ExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutpu
 
   for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
-    (*it)->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
+    (*it)->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 
   if (!alreadyWrittenAsTefTerm(symb_id, tef_terms))
     {
@@ -6367,7 +6520,7 @@ ExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutpu
             output << "TEF_" << indx << " = ";
 
           output << datatree.symbol_table.getName(symb_id) << "(";
-          writeExternalFunctionArguments(output, output_type, temporary_terms, tef_terms);
+          writeExternalFunctionArguments(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ");" << endl;
         }
     }
@@ -6523,6 +6676,7 @@ FirstDerivExternalFunctionNode::writeJsonOutput(ostream &output,
 void
 FirstDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                                             const temporary_terms_t &temporary_terms,
+                                            const temporary_terms_idxs_t &temporary_terms_idxs,
                                             deriv_node_temp_terms_t &tef_terms) const
 {
   assert(output_type != oMatlabOutsideModel);
@@ -6531,7 +6685,7 @@ FirstDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
     {
       output << "\\frac{\\partial " << datatree.symbol_table.getTeXName(symb_id)
              << "}{\\partial " << inputIndex << "}(";
-      writeExternalFunctionArguments(output, output_type, temporary_terms, tef_terms);
+      writeExternalFunctionArguments(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ")";
       return;
     }
@@ -6543,7 +6697,12 @@ FirstDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
       if (output_type == oMatlabDynamicModelSparse)
         output << "T" << idx << "(it_)";
       else
-        output << "T" << idx;
+        if (temporary_terms_idxs.empty() || IS_C(output_type))
+          output << "T" << idx;
+        else
+          output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                 << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type);
       return;
     }
 
@@ -6566,6 +6725,15 @@ FirstDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
            << LEFT_ARRAY_SUBSCRIPT(output_type) << tmpIndx << RIGHT_ARRAY_SUBSCRIPT(output_type);
 }
 
+void
+FirstDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                                            const temporary_terms_t &temporary_terms,
+                                            deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 void
 FirstDerivExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number,
                                         bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -6607,6 +6775,7 @@ FirstDerivExternalFunctionNode::compile(ostream &CompileCode, unsigned int &inst
 void
 FirstDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                                             const temporary_terms_t &temporary_terms,
+                                                            const temporary_terms_idxs_t &temporary_terms_idxs,
                                                             deriv_node_temp_terms_t &tef_terms) const
 {
   assert(output_type != oMatlabOutsideModel);
@@ -6618,7 +6787,7 @@ FirstDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Exp
   if (first_deriv_symb_id == symb_id)
     {
       expr_t parent = datatree.AddExternalFunction(symb_id, arguments);
-      parent->writeExternalFunctionOutput(output, output_type, temporary_terms,
+      parent->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs,
                                           tef_terms);
       return;
     }
@@ -6652,7 +6821,7 @@ FirstDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Exp
             output << "mxSetCell(prhs" << ending.str() << "[2], "
                    << i++ << ", "
                    << "mxCreateDoubleScalar("; // All external_function arguments are scalars
-            (*it)->writeOutput(output, output_type, temporary_terms, tef_terms);
+            (*it)->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
             output << "));" << endl;
           }
 
@@ -6699,7 +6868,7 @@ FirstDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Exp
                  << " = " << datatree.symbol_table.getName(first_deriv_symb_id) << "(";
         }
 
-      writeExternalFunctionArguments(output, output_type, temporary_terms, tef_terms);
+      writeExternalFunctionArguments(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 
       if (first_deriv_symb_id == eExtFunNotSet)
         output << "}";
@@ -6907,6 +7076,7 @@ SecondDerivExternalFunctionNode::writeJsonOutput(ostream &output,
 void
 SecondDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                                              const temporary_terms_t &temporary_terms,
+                                             const temporary_terms_idxs_t &temporary_terms_idxs,
                                              deriv_node_temp_terms_t &tef_terms) const
 {
   assert(output_type != oMatlabOutsideModel);
@@ -6915,7 +7085,7 @@ SecondDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
     {
       output << "\\frac{\\partial^2 " << datatree.symbol_table.getTeXName(symb_id)
              << "}{\\partial " << inputIndex1 << "\\partial " << inputIndex2 << "}(";
-      writeExternalFunctionArguments(output, output_type, temporary_terms, tef_terms);
+      writeExternalFunctionArguments(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
       output << ")";
       return;
     }
@@ -6927,7 +7097,12 @@ SecondDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
       if (output_type == oMatlabDynamicModelSparse)
         output << "T" << idx << "(it_)";
       else
-        output << "T" << idx;
+        if (temporary_terms_idxs.empty() || IS_C(output_type))
+          output << "T" << idx;
+        else
+          output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                 << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type);
       return;
     }
 
@@ -6962,9 +7137,19 @@ SecondDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType
              << LEFT_ARRAY_SUBSCRIPT(output_type) << tmpIndex1 << "," << tmpIndex2 << RIGHT_ARRAY_SUBSCRIPT(output_type);
 }
 
+void
+SecondDerivExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                                            const temporary_terms_t &temporary_terms,
+                                            deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 void
 SecondDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                                              const temporary_terms_t &temporary_terms,
+                                                             const temporary_terms_idxs_t &temporary_terms_idxs,
                                                              deriv_node_temp_terms_t &tef_terms) const
 {
   assert(output_type != oMatlabOutsideModel);
@@ -6976,7 +7161,7 @@ SecondDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Ex
   if (second_deriv_symb_id == symb_id)
     {
       expr_t parent = datatree.AddExternalFunction(symb_id, arguments);
-      parent->writeExternalFunctionOutput(output, output_type, temporary_terms,
+      parent->writeExternalFunctionOutput(output, output_type, temporary_terms, temporary_terms_idxs,
                                           tef_terms);
       return;
     }
@@ -7011,7 +7196,7 @@ SecondDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Ex
             output << "mxSetCell(prhs" << ending.str() << "[3], "
                    << i++ << ", "
                    << "mxCreateDoubleScalar("; // All external_function arguments are scalars
-            (*it)->writeOutput(output, output_type, temporary_terms, tef_terms);
+            (*it)->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
             output << "));" << endl;
           }
 
@@ -7060,7 +7245,7 @@ SecondDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Ex
                  << " = " << datatree.symbol_table.getName(second_deriv_symb_id) << "(";
         }
 
-      writeExternalFunctionArguments(output, output_type, temporary_terms, tef_terms);
+      writeExternalFunctionArguments(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
 
       if (second_deriv_symb_id == eExtFunNotSet)
         output << "}";
@@ -7218,6 +7403,7 @@ VarExpectationNode::cloneDynamic(DataTree &dynamic_datatree) const
 void
 VarExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                                 const temporary_terms_t &temporary_terms,
+                                const temporary_terms_idxs_t &temporary_terms_idxs,
                                 deriv_node_temp_terms_t &tef_terms) const
 {
   assert(output_type != oMatlabOutsideModel);
@@ -7237,13 +7423,27 @@ VarExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       if (output_type == oMatlabDynamicModelSparse)
         output << "T" << idx << "(it_)";
       else
-        output << "T" << idx;
+        if (temporary_terms_idxs.empty() || IS_C(output_type))
+          output << "T" << idx;
+        else
+          output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                 << temporary_terms_idxs[distance(temporary_terms.begin(), it)] + 1
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type);
       return;
     }
 
   output << "dynamic_var_forecast_" << model_name << "_" << forecast_horizon << "(" << yidx + 1 << ")";
 }
 
+void
+VarExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                                            const temporary_terms_t &temporary_terms,
+                                            deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 int
 VarExpectationNode::maxEndoLead() const
 {
@@ -7337,6 +7537,11 @@ VarExpectationNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, i
 {
 }
 
+void
+VarExpectationNode::collectModelLocalVariables(map<int, expr_t> &result) const
+{
+}
+
 void
 VarExpectationNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
 {
@@ -7580,6 +7785,7 @@ PacExpectationNode::cloneDynamic(DataTree &dynamic_datatree) const
 void
 PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                                 const temporary_terms_t &temporary_terms,
+                                const temporary_terms_idxs_t &temporary_terms_idxs,
                                 deriv_node_temp_terms_t &tef_terms) const
 {
   assert(output_type != oMatlabOutsideModel);
@@ -7673,6 +7879,15 @@ PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
   output << "];" << endl;
 }
 
+void
+PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                                            const temporary_terms_t &temporary_terms,
+                                            deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_idxs_t temporary_terms_idxs;
+  writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+}
+
 int
 PacExpectationNode::maxEndoLead() const
 {
@@ -7758,6 +7973,11 @@ PacExpectationNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, i
 {
 }
 
+void
+PacExpectationNode::collectModelLocalVariables(map<int, expr_t> &result) const
+{
+}
+
 void
 PacExpectationNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
 {
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 4632653508f6d0b05c0c62ddcb3fe5bb05c60c06..d3b2fdcbb13c6d8c8d77060872af70a8b74c8312 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -45,6 +45,8 @@ struct ExprNodeLess;
 //! Type for set of temporary terms
 /*! They are ordered by index number thanks to ExprNodeLess */
 typedef set<expr_t, ExprNodeLess> temporary_terms_t;
+/*! Keeps track of array indices of temporary_terms for writing */
+typedef vector<int> temporary_terms_idxs_t;
 
 //! set of temporary terms used in a block
 typedef set<int> temporary_terms_inuse_t;
@@ -219,6 +221,7 @@ class ExprNode
         \param[in,out] tef_terms the set of already written external function nodes
       */
       virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const = 0;
+      virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const = 0;
 
       //! returns true if the expr node contains an external function
       virtual bool containsExternalFunction() const = 0;
@@ -231,6 +234,7 @@ class ExprNode
 
       //! Writes output of node, using a Txxx notation for nodes in temporary_terms
       void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
+      void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs) const;
 
       //! Writes output of node in JSON syntax
       virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic = true) const = 0;
@@ -240,6 +244,7 @@ class ExprNode
       //! Writes the output for an external function, ensuring that the external function is called as few times as possible using temporary terms
       virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                                const temporary_terms_t &temporary_terms,
+                                               const temporary_terms_idxs_t &temporary_terms_idxs,
                                                deriv_node_temp_terms_t &tef_terms) const;
 
       //! Write the JSON output of an external function in a string vector
@@ -263,6 +268,9 @@ class ExprNode
       */
       virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const = 0;
 
+      //! Used in the collection of Model Local Variables for writing as Temporary Terms in Static and Dynamic files
+      virtual void collectModelLocalVariables(map<int, expr_t> &result) const = 0;
+
       //! Computes the set of all variables of a given symbol type in the expression (without information on lags)
       /*!
         Variables are stored as symb_id.
@@ -539,9 +547,11 @@ public:
   };
   virtual void prepareForDerivation();
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
+  virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
   virtual bool containsExternalFunction() const;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual void collectModelLocalVariables(map<int, expr_t> &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, 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;
@@ -601,9 +611,11 @@ public:
   VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg);
   virtual void prepareForDerivation();
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
+  virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
   virtual bool containsExternalFunction() const;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual void collectModelLocalVariables(map<int, expr_t> &result) const;
   virtual void computeTemporaryTerms(map<expr_t, int > &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      map<expr_t, pair<int, int> > &first_occurence,
@@ -695,10 +707,12 @@ public:
                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) 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 writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
   virtual bool containsExternalFunction() const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs,
                                            deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonExternalFunctionOutput(vector<string> &efout,
                                                const temporary_terms_t &temporary_terms,
@@ -715,6 +729,7 @@ public:
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual void collectModelLocalVariables(map<int, expr_t> &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, EvalExternalFunctionException);
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
@@ -801,10 +816,12 @@ public:
                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) 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 writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
   virtual bool containsExternalFunction() const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs,
                                            deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonExternalFunctionOutput(vector<string> &efout,
                                                const temporary_terms_t &temporary_terms,
@@ -821,6 +838,7 @@ public:
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual void collectModelLocalVariables(map<int, expr_t> &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, int derivOrder) throw (EvalException, EvalExternalFunctionException);
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
@@ -926,10 +944,12 @@ public:
                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) 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 writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
   virtual bool containsExternalFunction() const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs,
                                            deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonExternalFunctionOutput(vector<string> &efout,
                                                const temporary_terms_t &temporary_terms,
@@ -946,6 +966,7 @@ public:
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual void collectModelLocalVariables(map<int, expr_t> &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, EvalExternalFunctionException);
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
@@ -1012,7 +1033,7 @@ protected:
   //! Returns the index in the tef_terms map of this external function
   int getIndxInTefTerms(int the_symb_id, deriv_node_temp_terms_t &tef_terms) const throw (UnknownFunctionNameAndArgs);
   //! Helper function to write output arguments of any given external function
-  void writeExternalFunctionArguments(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
+  void writeExternalFunctionArguments(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   void writeJsonExternalFunctionArguments(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
 public:
   AbstractExternalFunctionNode(DataTree &datatree_arg, int symb_id_arg,
@@ -1022,10 +1043,12 @@ public:
                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) const = 0;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const = 0;
+  virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const = 0;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic = true) const = 0;
   virtual bool containsExternalFunction() const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs,
                                            deriv_node_temp_terms_t &tef_terms) const = 0;
   virtual void writeJsonExternalFunctionOutput(vector<string> &efout,
                                                const temporary_terms_t &temporary_terms,
@@ -1042,6 +1065,7 @@ public:
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const = 0;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual void collectModelLocalVariables(map<int, expr_t> &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, EvalExternalFunctionException);
   unsigned int compileExternalFunctionArguments(ostream &CompileCode, unsigned int &instruction_number,
@@ -1105,9 +1129,11 @@ public:
                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) 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 writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs,
                                            deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonExternalFunctionOutput(vector<string> &efout,
                                                const temporary_terms_t &temporary_terms,
@@ -1150,6 +1176,7 @@ 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 writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -1157,6 +1184,7 @@ public:
                        deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs,
                                            deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonExternalFunctionOutput(vector<string> &efout,
                                                const temporary_terms_t &temporary_terms,
@@ -1194,6 +1222,7 @@ 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 writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -1201,6 +1230,7 @@ public:
                        deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
+                                           const temporary_terms_idxs_t &temporary_terms_idxs,
                                            deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeJsonExternalFunctionOutput(vector<string> &efout,
                                                const temporary_terms_t &temporary_terms,
@@ -1229,6 +1259,7 @@ public:
                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) 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 writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, 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,
@@ -1266,6 +1297,7 @@ public:
                        deriv_node_temp_terms_t &tef_terms) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual void collectModelLocalVariables(map<int, expr_t> &result) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
   virtual bool isDiffPresent(void) const;
@@ -1306,6 +1338,7 @@ public:
                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) 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 writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, 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,
@@ -1343,6 +1376,7 @@ public:
                        deriv_node_temp_terms_t &tef_terms) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual void collectModelLocalVariables(map<int, expr_t> &result) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
   virtual bool isDiffPresent(void) const;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index aea4f6b70b824b38563e37f82ebcfc74b267613a..e9ea97e4c8f59dcce646c230a28d6b488934c184 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -1120,10 +1120,25 @@ ModelTree::computeTemporaryTerms(bool is_matlab)
 {
   map<expr_t, pair<int, NodeTreeReference> > reference_count;
   temporary_terms.clear();
+  temporary_terms_mlv.clear();
   temporary_terms_res.clear();
   temporary_terms_g1.clear();
   temporary_terms_g2.clear();
   temporary_terms_g3.clear();
+
+  // Collect all model local variables appearing in equations. See #101
+  // All used model local variables are automatically set as temporary variables
+  map<int, expr_t> used_local_vars;
+  for (size_t i = 0; i < equations.size(); i++)
+    equations[i]->collectModelLocalVariables(used_local_vars);
+
+  for (map<int, expr_t>::const_iterator it = used_local_vars.begin();
+       it != used_local_vars.end(); it++)
+    {
+      temporary_terms_mlv[it->second] = local_variables_table.find(it->first)->second;
+      reference_count[it->second] = make_pair(MIN_COST(is_matlab)+1, eResiduals);
+    }
+
   map<NodeTreeReference, temporary_terms_t> temp_terms_map;
   temp_terms_map[eResiduals] = temporary_terms_res;
   temp_terms_map[eFirstDeriv] = temporary_terms_g1;
@@ -1154,20 +1169,54 @@ ModelTree::computeTemporaryTerms(bool is_matlab)
                                       temp_terms_map,
                                       is_matlab, eThirdDeriv);
 
+  int idx = 0;
   for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.begin();
        it != temp_terms_map.end(); it++)
-    temporary_terms.insert(it->second.begin(), it->second.end());
+    {
+      temporary_terms.insert(it->second.begin(), it->second.end());
+      temporary_terms_idxs.push_back(idx++);
+    }
 
   temporary_terms_res = temp_terms_map[eResiduals];
   temporary_terms_g1  = temp_terms_map[eFirstDeriv];
   temporary_terms_g2  = temp_terms_map[eSecondDeriv];
   temporary_terms_g3  = temp_terms_map[eThirdDeriv];
+
+  temporary_terms_mlv_idxs.clear();
+  temporary_terms_res_idxs.clear();
+  temporary_terms_g1_idxs.clear();
+  temporary_terms_g2_idxs.clear();
+  temporary_terms_g3_idxs.clear();
+
+  idx = 0;
+  for (map<expr_t, expr_t>::const_iterator it = temporary_terms_mlv.begin();
+       it != temporary_terms_mlv.end(); it++)
+    temporary_terms_mlv_idxs.push_back(idx++);
+
+  for (temporary_terms_t::const_iterator it = temporary_terms_res.begin();
+       it != temporary_terms_res.end(); it++)
+    temporary_terms_res_idxs.push_back(idx++);
+
+  for (temporary_terms_t::const_iterator it = temporary_terms_g1.begin();
+       it != temporary_terms_g1.end(); it++)
+    temporary_terms_g1_idxs.push_back(idx++);
+
+  for (temporary_terms_t::const_iterator it = temporary_terms_g2.begin();
+       it != temporary_terms_g2.end(); it++)
+    temporary_terms_g2_idxs.push_back(idx++);
+
+  for (temporary_terms_t::const_iterator it = temporary_terms_g3.begin();
+       it != temporary_terms_g3.end(); it++)
+    temporary_terms_g3_idxs.push_back(idx++);
 }
 
 void
 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
 {
+  temporary_terms_idxs_t tti, ttm1i;
+  writeTemporaryTerms(tt, tti, ttm1, ttm1i, output, output_type, tef_terms);
+  /*
   // Local var used to keep track of temp nodes already written
   temporary_terms_t tt2 = ttm1;
   for (temporary_terms_t::const_iterator it = tt.begin();
@@ -1193,6 +1242,72 @@ ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, const temporary_term
         // Insert current node into tt2
         tt2.insert(*it);
       }
+  */
+}
+
+void
+ModelTree::writeModelLocalVariableTemporaryTerms(const temporary_terms_t &tto, const map<expr_t, expr_t> &tt,
+                                                 const temporary_terms_idxs_t &ttidxs,
+                                                 ostream &output, ExprNodeOutputType output_type,
+                                                 deriv_node_temp_terms_t &tef_terms) const
+{
+  temporary_terms_t tt2;
+  temporary_terms_idxs_t tt2idxs;
+  int ttidx = 0;
+  for (map<expr_t, expr_t>::const_iterator it = tt.begin();
+       it != tt.end(); it++, ttidx++)
+    {
+      if (IS_C(output_type))
+        output << "double ";
+      else if (IS_JULIA(output_type))
+        output << "    @inbounds const ";
+
+      it->first->writeOutput(output, output_type, tto, ttidxs, tef_terms);
+      output << " = ";
+      it->second->writeOutput(output, output_type, tt2, tt2idxs, tef_terms);
+
+      if (IS_C(output_type) || IS_MATLAB(output_type))
+        output << ";";
+      output << endl;
+
+      // Insert current node into tt2
+      tt2.insert(it->first);
+      tt2idxs.push_back(ttidxs[ttidx]);
+    }
+}
+
+void
+ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, const temporary_terms_idxs_t &ttidxs,
+                               const temporary_terms_t &ttm1, const temporary_terms_idxs_t &ttm1idxs,
+                               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 = ttm1;
+  temporary_terms_idxs_t tt2idxs = ttm1idxs;
+  int ttidx = 0;
+  for (temporary_terms_t::const_iterator it = tt.begin();
+       it != tt.end(); it++, ttidx++)
+    {
+      if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
+        (*it)->writeExternalFunctionOutput(output, output_type, tt2, tt2idxs, tef_terms);
+
+      if (IS_C(output_type))
+        output << "double ";
+      else if (IS_JULIA(output_type))
+        output << "    @inbounds ";
+
+      (*it)->writeOutput(output, output_type, tt, ttidxs, tef_terms);
+      output << " = ";
+      (*it)->writeOutput(output, output_type, tt2, tt2idxs, tef_terms);
+
+      if (IS_C(output_type) || IS_MATLAB(output_type))
+        output << ";";
+      output << endl;
+
+      // Insert current node into tt2
+      tt2.insert(*it);
+      tt2idxs.push_back(ttidxs[ttidx]);
+    }
 }
 
 void
@@ -1407,14 +1522,13 @@ ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_n
 void
 ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const
 {
+  return;
   /* Collect all model local variables appearing in equations, and print only
      them. Printing unused model local variables can lead to a crash (see
      ticket #101). */
   set<int> used_local_vars;
-
-  // Use an empty set for the temporary terms
   const temporary_terms_t tt;
-
+  const temporary_terms_idxs_t tti;
   for (size_t i = 0; i < equations.size(); i++)
     equations[i]->collectVariables(eModelLocalVariable, used_local_vars);
 
@@ -1424,17 +1538,17 @@ ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_t
       {
         int id = *it;
         expr_t value = local_variables_table.find(id)->second;
-        value->writeExternalFunctionOutput(output, output_type, tt, tef_terms);
+        value->writeExternalFunctionOutput(output, output_type, tt, tti, tef_terms);
 
         if (IS_C(output_type))
           output << "double ";
         else if (IS_JULIA(output_type))
-          output << "  @inbounds ";
+          output << "    @inbounds ";
 
         /* We append underscores to avoid name clashes with "g1" or "oo_" (see
            also VariableNode::writeOutput) */
         output << symbol_table.getName(id) << "__ = ";
-        value->writeOutput(output, output_type, tt, tef_terms);
+        value->writeOutput(output, output_type, tt, tti, tef_terms);
         output << ";" << endl;
       }
 }
@@ -1491,12 +1605,16 @@ ModelTree::writeJsonModelLocalVariables(ostream &output, deriv_node_temp_terms_t
 void
 ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const
 {
-  temporary_terms_t temp_terms;
-  if (IS_JULIA(output_type))
-    temp_terms = temporary_terms_res;
-  else
-    temp_terms = temporary_terms;
+  temporary_terms_t tt;
+  temporary_terms_idxs_t ttidxs;
+  writeModelEquations(output, output_type, tt, ttidxs);
+}
 
+void
+ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type,
+                               const temporary_terms_t &temporary_terms,
+                               const temporary_terms_idxs_t &temporary_terms_idxs) const
+{
   for (int eq = 0; eq < (int) equations.size(); eq++)
     {
       BinaryOpNode *eq_node = equations[eq];
@@ -1514,35 +1632,39 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type)
         }
 
       if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
-        {
-          if (IS_JULIA(output_type))
-            output << "  @inbounds ";
-          output << "lhs =";
-          lhs->writeOutput(output, output_type, temp_terms);
-          output << ";" << endl;
-
-          if (IS_JULIA(output_type))
-            output << "  @inbounds ";
-          output << "rhs =";
-          rhs->writeOutput(output, output_type, temp_terms);
-          output << ";" << endl;
-
-          if (IS_JULIA(output_type))
-            output << "  @inbounds ";
-          output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
-                 << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
-                 << RIGHT_ARRAY_SUBSCRIPT(output_type)
-                 << "= lhs-rhs;" << endl;
-        }
+        if (IS_JULIA(output_type))
+          {
+            output << "    @inbounds residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                   << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                   << " = (";
+            lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+            output << ") - (";
+            rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+            output << ")" << endl;
+          }
+        else
+          {
+            output << "lhs = ";
+            lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+            output << ";" << endl
+                   << "rhs = ";
+            rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+            output << ";" << endl
+                   << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                   << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type)
+                   << " = lhs - rhs;" << endl;
+          }
       else // The right hand side of the equation is empty ==> residual=lhs;
         {
           if (IS_JULIA(output_type))
-            output << "  @inbounds ";
+            output << "    @inbounds ";
           output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
                  << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
                  << RIGHT_ARRAY_SUBSCRIPT(output_type)
                  << " = ";
-          lhs->writeOutput(output, output_type, temp_terms);
+          lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
           output << ";" << endl;
         }
     }
@@ -1787,9 +1909,8 @@ ModelTree::set_cutoff_to_zero()
 void
 ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
 {
-  output << "  ";
   if (IS_JULIA(output_type))
-    output << "@inbounds ";
+    output << "    @inbounds ";
   output << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
   if (IS_MATLAB(output_type) || IS_JULIA(output_type))
     output << eq_nb + 1 << "," << col_nb + 1;
@@ -1801,7 +1922,7 @@ ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutput
 void
 ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
 {
-  output << "  v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
+  output << "v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
   if (IS_MATLAB(output_type) || IS_JULIA(output_type))
     output << row_nb + 1 << "," << col_nb + 1;
   else
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 312169a1f90bfc7d3af1d8d3f2f14f76cfe6f1ed..a0553a004401051184907ed40f8b37b805659b2f 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -134,13 +134,21 @@ protected:
   */
   third_derivatives_t hessian_params_derivatives;
 
-  //! Temporary terms for the static/dynamic file (those which will be noted Txxxx)
+  //! Temporary terms for the static/dynamic file (those which will be noted T[x])
   temporary_terms_t temporary_terms;
+  map<expr_t, expr_t> temporary_terms_mlv;
   temporary_terms_t temporary_terms_res;
   temporary_terms_t temporary_terms_g1;
   temporary_terms_t temporary_terms_g2;
   temporary_terms_t temporary_terms_g3;
 
+  temporary_terms_idxs_t temporary_terms_idxs;
+  temporary_terms_idxs_t temporary_terms_mlv_idxs;
+  temporary_terms_idxs_t temporary_terms_res_idxs;
+  temporary_terms_idxs_t temporary_terms_g1_idxs;
+  temporary_terms_idxs_t temporary_terms_g2_idxs;
+  temporary_terms_idxs_t temporary_terms_g3_idxs;
+
   //! Temporary terms for the file containing parameters derivatives
   temporary_terms_t params_derivs_temporary_terms;
   temporary_terms_t params_derivs_temporary_terms_res;
@@ -183,6 +191,9 @@ protected:
   void computeParamsDerivativesTemporaryTerms();
   //! Writes temporary terms
   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;
+  void writeTemporaryTerms(const temporary_terms_t &tt, const temporary_terms_idxs_t &ttidxs,
+                           const temporary_terms_t &ttm1, const temporary_terms_idxs_t &ttm1idxs,
+                           ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
   void writeJsonTemporaryTerms(const temporary_terms_t &tt, const temporary_terms_t &ttm1, ostream &output, deriv_node_temp_terms_t &tef_terms, string &concat) 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;
@@ -195,8 +206,15 @@ protected:
   //! Writes model local variables
   /*! No temporary term is used in the output, so that local parameters declarations can be safely put before temporary terms declaration in the output files */
   void writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
+  void writeModelLocalVariableTemporaryTerms(const temporary_terms_t &tto, const map<expr_t, expr_t> &tt,
+                                             const temporary_terms_idxs_t &ttidxs,
+                                             ostream &output, ExprNodeOutputType output_type,
+                                             deriv_node_temp_terms_t &tef_terms) const;
   //! Writes model equations
   void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const;
+  void writeModelEquations(ostream &output, ExprNodeOutputType output_type,
+                           const temporary_terms_t &temporary_terms,
+                           const temporary_terms_idxs_t &temporary_terms_idxs) const;
   //! Writes JSON model equations
   //! if residuals = true, we are writing the dynamic/static model.
   //! Otherwise, just the model equations (with line numbers, no tmp terms)
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 2391de980797505f7564486e5472c2bf258d35d9..bd67ebd28e115a514d06f78f28f1168ceb897eb5 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -274,6 +274,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
         output << "  residual=zeros(" << block_mfs << ",1);\n";
 
       // The equations
+      temporary_terms_idxs_t temporary_terms_idxs;
       for (unsigned int i = 0; i < block_size; i++)
         {
           if (!global_temporary_terms)
@@ -287,12 +288,12 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
                    it != v_temporary_terms[block][i].end(); it++)
                 {
                   if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
-                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, tef_terms);
+                    (*it)->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
 
                   output << "  " <<  sps;
-                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, tef_terms);
+                  (*it)->writeOutput(output, local_output_type, local_temporary_terms, temporary_terms_idxs, tef_terms);
                   output << " = ";
-                  (*it)->writeOutput(output, local_output_type, tt2, tef_terms);
+                  (*it)->writeOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
                   // Insert current node into tt2
                   tt2.insert(*it);
                   output << ";" << endl;
@@ -1161,84 +1162,202 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
 void
 StaticModel::writeStaticMFile(const string &func_name) const
 {
-  // Writing comments and function definition command
-  string filename = func_name + "_static.m";
+  writeStaticModel(func_name, false, false);
+}
 
+void
+StaticModel::writeWrapperFunctions(const string &basename, const string &ending) const
+{
+  string name;
+  if (ending == "g1")
+    name = basename + "_resid_g1";
+  else if (ending == "g2")
+    name= basename + "_resid_g1_g2";
+  else if (ending == "g3")
+    name = basename + "_resid_g1_g2_g3";
+
+  string filename = name + ".m";
   ofstream output;
   output.open(filename.c_str(), ios::out | ios::binary);
   if (!output.is_open())
     {
-      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  if (ending == "g1")
+    output << "function [residual, g1] = " << name << "(T, y, x, params)" << endl
+           << "% function [residual, g1] = " << name << "(T, y, x, params)" << endl;
+  else if (ending == "g2")
+    output << "function [residual, g1, g2] = " << name << "(T, y, x, params)" << endl
+           << "% function [residual, g1, g2] = " << name << "(T, y, x, params)" << endl;
+  else if (ending == "g3")
+    output << "function [residual, g1, g2, g3] = " << name << "(T, y, x, params)" << endl
+           << "% function [residual, g1, g2, g3] = " << name << "(T, y, x, params)" << endl;
+
+  output << "%" << endl
+         << "% Wrapper function automatically created by Dynare" << endl
+         << "%" << endl << endl;
+
+  if (ending == "g1")
+    output << "    T = " << basename + "_" + ending + "_tt(T, y, x, params);" << endl
+           << "    residual = " << basename + "_resid(T, y, x, params, false);" << endl
+           << "    g1       = " << basename + "_g1(T, y, x, params, false);" << endl;
+  else if (ending == "g2")
+    output << "    T = " << basename + "_" + ending + "_tt(T, y, x, params);" << endl
+           << "    [residual, g1] = " << basename + "_resid_g1(T, y, x, params, false);" << endl
+           << "    g2       = " << basename + "_g2(T, y, x, params, false);" << endl;
+  else if (ending == "g3")
+    output << "    T = " << basename + "_" + ending + "_tt(T, y, x, params);" << endl
+           << "    [residual, g1, g2] = " << basename + "_resid_g1(T, y, x, params, false);" << endl
+           << "    g3       = " << basename + "_g3(T, y, x, params, false);" << endl;
+
+  output << endl << "end" << endl;
+  output.close();
+}
+
+void
+StaticModel::writeStaticModelHelper(const string &name, const string &retvalname,
+                                    const string &name_tt, size_t ttlen,
+                                    const string &previous_tt_name,
+                                    const ostringstream &init_s, const ostringstream &end_s,
+                                    const ostringstream &s, const ostringstream &s_tt) const
+{
+  string filename =  name_tt + ".m";
+  ofstream output;
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
 
-  output << "function [residual, g1, g2, g3] = " << func_name + "_static(y, x, params)" << endl
+  output << "function T = " << name_tt << "(T, y, x, params)" << endl
+         << "% function T = " << name_tt << "(T, y, x, params)" << endl
+         << "%" << endl
+         << "% File created by Dynare Preprocessor from .mod file" << endl
          << "%" << endl
-         << "% Status : Computes static model for Dynare" << endl
+         << "% Inputs:" << endl
+         << "%   T         [#temp variables by 1]  double   vector of temporary terms to be filled by function" << endl
+         << "%   y         [M_.endo_nbr by 1]      double   vector of endogenous variables in declaration order" << endl
+         << "%   x         [M_.exo_nbr by 1]       double   vector of exogenous variables in declaration order" << endl
+         << "%   params    [M_.param_nbr by 1]     double   vector of parameter values in declaration order" << endl
          << "%" << endl
-         << "% Inputs : " << endl
-         << "%   y         [M_.endo_nbr by 1] double    vector of endogenous variables in declaration order" << endl
-         << "%   x         [M_.exo_nbr by 1] double     vector of exogenous variables in declaration order" << endl
-         << "%   params    [M_.param_nbr by 1] double   vector of parameter values in declaration order" << endl
+         << "% Output:" << endl
+         << "%   T         [#temp variables by 1]  double   vector of temporary terms" << endl
+         << "%" << endl << endl
+         << "assert(length(T) >= " << ttlen << ");" << endl
+         << endl;
+
+  if (!previous_tt_name.empty())
+    output << "T = " << previous_tt_name << "(T, y, x, params);" << endl << endl;
+
+  output << s_tt.str() << endl
+         << "end" << endl;
+  output.close();
+
+  filename = name + ".m";
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  output << "function " << retvalname << " = " << name << "(T, y, x, params, T_flag)" << endl
+         << "% function " << retvalname << " = " << name << "(T, y, x, params, T_flag)" << endl
          << "%" << endl
-         << "% Outputs:" << endl
-         << "%   residual  [M_.endo_nbr by 1] double    vector of residuals of the static model equations " << endl
-         << "%                                          in order of declaration of the equations." << endl
-         << "%                                          Dynare may prepend or append auxiliary equations, see M_.aux_vars" << endl
-         << "%   g1        [M_.endo_nbr by M_.endo_nbr] double    Jacobian matrix of the static model equations;" << endl
-         << "%                                                       columns: variables in declaration order" << endl
-         << "%                                                       rows: equations in order of declaration" << endl
-         << "%   g2        [M_.endo_nbr by (M_.endo_nbr)^2] double   Hessian matrix of the static model equations;" << endl
-         << "%                                                       columns: variables in declaration order" << endl
-         << "%                                                       rows: equations in order of declaration" << endl
-         << "%   g3        [M_.endo_nbr by (M_.endo_nbr)^3] double   Third derivatives matrix of the static model equations;" << endl
-         << "%                                                       columns: variables in declaration order" << endl
-         << "%                                                       rows: equations in order of declaration" << endl
+         << "% File created by Dynare Preprocessor from .mod file" << endl
          << "%" << endl
+         << "% Inputs:" << endl
+         << "%   T         [#temp variables by 1]  double   vector of temporary terms to be filled by function" << endl
+         << "%   y         [M_.endo_nbr by 1]      double   vector of endogenous variables in declaration order" << endl
+         << "%   x         [M_.exo_nbr by 1]       double   vector of exogenous variables in declaration order" << endl
+         << "%   params    [M_.param_nbr by 1]     double   vector of parameter values in declaration order" << endl
+         << "%                                              to evaluate the model" << endl
+         << "%   T_flag    boolean                 boolean  flag saying whether or not to calculate temporary terms" << endl
          << "%" << endl
-         << "% Warning : this file is generated automatically by Dynare" << endl
-         << "%           from model file (.mod)" << endl << endl;
-
-  writeStaticModel(output, false, false);
-  output << "end" << endl;
+         << "% Output:" << endl
+         << "%   " << retvalname << endl
+         << "%" << endl << endl;
+
+  if (!name_tt.empty())
+    output << "if T_flag" << endl
+           << "    T = " << name_tt << "(T, y, x, params);" << endl
+           << "end" << endl;
+
+  output << init_s.str() << endl
+         << s.str()
+         << end_s.str() << endl
+         << "end" << endl;
   output.close();
 }
 
 void
 StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const
 {
-  ostringstream model_local_vars_output;   // Used for storing model local vars
-  ostringstream model_output;              // Used for storing model
-  ostringstream jacobian_output;           // Used for storing jacobian equations
-  ostringstream hessian_output;            // Used for storing Hessian equations
-  ostringstream third_derivatives_output;  // Used for storing third order derivatives equations
+  writeStaticModel("", StaticOutput, use_dll, julia);
+}
+
+void
+StaticModel::writeStaticModel(const string &basename, bool use_dll, bool julia) const
+{
+  ofstream StaticOutput;
+  writeStaticModel(basename, StaticOutput, use_dll, julia);
+}
+
+void
+StaticModel::writeStaticModel(const string &basename,
+                              ostream &StaticOutput, bool use_dll, bool julia) const
+{
+  ostringstream model_tt_output;             // Used for storing model temp vars
+  ostringstream model_output;                // Used for storing model equations
+  ostringstream jacobian_tt_output;          // Used for storing jacobian temp vars
+  ostringstream jacobian_output;             // Used for storing jacobian equations
+  ostringstream hessian_tt_output;           // Used for storing Hessian temp vars
+  ostringstream hessian_output;              // Used for storing Hessian equations
+  ostringstream third_derivatives_tt_output; // Used for storing third order derivatives temp terms
+  ostringstream third_derivatives_output;    // Used for storing third order derivatives equations
   ostringstream for_sym;
+
   ExprNodeOutputType output_type = (use_dll ? oCStaticModel :
                                     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, temp_term_union_m_1, model_output, output_type, tef_terms);
-
-  writeModelEquations(model_output, output_type);
+  temporary_terms_t temp_term_union;
+  temporary_terms_idxs_t temp_term_union_idxs;
+
+  for (map<expr_t, expr_t>::const_iterator it = temporary_terms_mlv.begin();
+       it != temporary_terms_mlv.end(); it++)
+    temp_term_union.insert(it->first);
+  writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_mlv,
+                                        temporary_terms_mlv_idxs,
+                                        model_tt_output, output_type, tef_terms);
+  temp_term_union_idxs = temporary_terms_mlv_idxs;
+
+  writeTemporaryTerms(temporary_terms_res, temporary_terms_res_idxs,
+                      temp_term_union, temp_term_union_idxs,
+                      model_tt_output, output_type, tef_terms);
+  temp_term_union.insert(temporary_terms_res.begin(), temporary_terms_res.end());
+  temp_term_union_idxs.insert(temp_term_union_idxs.end(),
+                              temporary_terms_res_idxs.begin(), temporary_terms_res_idxs.end());
+
+  writeModelEquations(model_output, output_type, temp_term_union, temp_term_union_idxs);
 
   int nrows = equations.size();
   int JacobianColsNbr = symbol_table.endo_nbr();
   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, temp_term_empty, jacobian_output, output_type, tef_terms);
-    else
-      writeTemporaryTerms(temp_term_union, temp_term_union_m_1, jacobian_output, output_type, tef_terms);
+    {
+      writeTemporaryTerms(temporary_terms_g1, temporary_terms_g1_idxs,
+                          temp_term_union, temp_term_union_idxs,
+                          jacobian_tt_output, output_type, tef_terms);
+      temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end());
+      temp_term_union_idxs.insert(temp_term_union_idxs.end(),
+                                  temporary_terms_g1_idxs.begin(), temporary_terms_g1_idxs.end());
+    }
   for (first_derivatives_t::const_iterator it = first_derivatives.begin();
        it != first_derivatives.end(); it++)
     {
@@ -1248,156 +1367,163 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
 
       jacobianHelper(jacobian_output, eq, symbol_table.getTypeSpecificID(symb_id), output_type);
       jacobian_output << "=";
-      d1->writeOutput(jacobian_output, output_type, temp_term_union, tef_terms);
+      d1->writeOutput(jacobian_output, output_type,
+                      temp_term_union, temp_term_union_idxs, tef_terms);
       jacobian_output << ";" << endl;
     }
 
   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, temp_term_empty, hessian_output, output_type, tef_terms);
-    else
-      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++)
     {
-      int eq = it->first.first;
-      int symb_id1 = getSymbIDByDerivID(it->first.second.first);
-      int symb_id2 = getSymbIDByDerivID(it->first.second.second);
-      expr_t d2 = it->second;
+      writeTemporaryTerms(temporary_terms_g2, temporary_terms_g2_idxs,
+                          temp_term_union, temp_term_union_idxs,
+                          hessian_tt_output, output_type, tef_terms);
+      temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end());
+      temp_term_union_idxs.insert(temp_term_union_idxs.end(),
+                                  temporary_terms_g2_idxs.begin(), temporary_terms_g2_idxs.end());
 
-      int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
-      int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
+      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++)
+        {
+          int eq = it->first.first;
+          int symb_id1 = getSymbIDByDerivID(it->first.second.first);
+          int symb_id2 = getSymbIDByDerivID(it->first.second.second);
+          expr_t d2 = it->second;
 
-      int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
-      int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
+          int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
+          int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
 
-      if (output_type == oJuliaStaticModel)
-        {
-          for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
-          hessian_output << "  @inbounds " << for_sym.str() << " = ";
-          d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
-          hessian_output << endl;
-        }
-      else
-        {
-          sparseHelper(2, hessian_output, k, 0, output_type);
-          hessian_output << "=" << eq + 1 << ";" << endl;
+          int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
+          int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
 
-          sparseHelper(2, hessian_output, k, 1, output_type);
-          hessian_output << "=" << col_nb + 1 << ";" << endl;
+          if (output_type == oJuliaStaticModel)
+            {
+              for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
+              hessian_output << "  @inbounds " << for_sym.str() << " = ";
+              d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
+              hessian_output << endl;
+            }
+          else
+            {
+              sparseHelper(2, hessian_output, k, 0, output_type);
+              hessian_output << "=" << eq + 1 << ";" << endl;
 
-          sparseHelper(2, hessian_output, k, 2, output_type);
-          hessian_output << "=";
-          d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
-          hessian_output << ";" << endl;
+              sparseHelper(2, hessian_output, k, 1, output_type);
+              hessian_output << "=" << col_nb + 1 << ";" << endl;
 
-          k++;
-        }
+              sparseHelper(2, hessian_output, k, 2, output_type);
+              hessian_output << "=";
+              d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
+              hessian_output << ";" << endl;
 
-      // Treating symetric elements
-      if (symb_id1 != symb_id2)
-        if (output_type == oJuliaStaticModel)
-          hessian_output << "  @inbounds g2[" << eq + 1 << "," << col_nb_sym + 1 << "] = "
-                         << for_sym.str() << endl;
-        else
-          {
-            sparseHelper(2, hessian_output, k, 0, output_type);
-            hessian_output << "=" << eq + 1 << ";" << endl;
+              k++;
+            }
 
-            sparseHelper(2, hessian_output, k, 1, output_type);
-            hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
+          // Treating symetric elements
+          if (symb_id1 != symb_id2)
+            if (output_type == oJuliaStaticModel)
+              hessian_output << "  @inbounds g2[" << eq + 1 << "," << col_nb_sym + 1 << "] = "
+                             << for_sym.str() << endl;
+            else
+              {
+                sparseHelper(2, hessian_output, k, 0, output_type);
+                hessian_output << "=" << eq + 1 << ";" << endl;
 
-            sparseHelper(2, hessian_output, k, 2, output_type);
-            hessian_output << "=";
-            sparseHelper(2, hessian_output, k-1, 2, output_type);
-            hessian_output << ";" << endl;
+                sparseHelper(2, hessian_output, k, 1, output_type);
+                hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
 
-            k++;
-          }
+                sparseHelper(2, hessian_output, k, 2, output_type);
+                hessian_output << "=";
+                sparseHelper(2, hessian_output, k-1, 2, output_type);
+                hessian_output << ";" << endl;
+
+                k++;
+              }
+        }
     }
 
   // 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, temp_term_empty, third_derivatives_output, output_type, tef_terms);
-    else
-      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++)
     {
-      int eq = it->first.first;
-      int var1 = it->first.second.first;
-      int var2 = it->first.second.second.first;
-      int var3 = it->first.second.second.second;
-      expr_t d3 = it->second;
-
-      int id1 = getSymbIDByDerivID(var1);
-      int id2 = getSymbIDByDerivID(var2);
-      int id3 = getSymbIDByDerivID(var3);
-
-      // Reference column number for the g3 matrix
-      int ref_col = id1 * hessianColsNbr + id2 * JacobianColsNbr + id3;
+      writeTemporaryTerms(temporary_terms_g3, temporary_terms_g3_idxs,
+                          temp_term_union, temp_term_union_idxs,
+                          third_derivatives_tt_output, output_type, tef_terms);
+      temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end());
+      temp_term_union_idxs.insert(temp_term_union_idxs.end(),
+                                  temporary_terms_g3_idxs.begin(), temporary_terms_g3_idxs.end());
 
-      if (output_type == oJuliaStaticModel)
+      int 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++)
         {
-          for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
-          third_derivatives_output << "  @inbounds " << for_sym.str() << " = ";
-          d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
-          third_derivatives_output << endl;
-        }
-      else
-        {
-          sparseHelper(3, third_derivatives_output, k, 0, output_type);
-          third_derivatives_output << "=" << eq + 1 << ";" << endl;
-
-          sparseHelper(3, third_derivatives_output, k, 1, output_type);
-          third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
+          int eq = it->first.first;
+          int var1 = it->first.second.first;
+          int var2 = it->first.second.second.first;
+          int var3 = it->first.second.second.second;
+          expr_t d3 = it->second;
 
-          sparseHelper(3, third_derivatives_output, k, 2, output_type);
-          third_derivatives_output << "=";
-          d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
-          third_derivatives_output << ";" << endl;
-        }
+          int id1 = getSymbIDByDerivID(var1);
+          int id2 = getSymbIDByDerivID(var2);
+          int id3 = getSymbIDByDerivID(var3);
 
-      // Compute the column numbers for the 5 other permutations of (id1,id2,id3)
-      // and store them in a set (to avoid duplicates if two indexes are equal)
-      set<int> cols;
-      cols.insert(id1 * hessianColsNbr + id3 * JacobianColsNbr + id2);
-      cols.insert(id2 * hessianColsNbr + id1 * JacobianColsNbr + id3);
-      cols.insert(id2 * hessianColsNbr + id3 * JacobianColsNbr + id1);
-      cols.insert(id3 * hessianColsNbr + id1 * JacobianColsNbr + id2);
-      cols.insert(id3 * hessianColsNbr + id2 * JacobianColsNbr + id1);
+          // Reference column number for the g3 matrix
+          int ref_col = id1 * hessianColsNbr + id2 * JacobianColsNbr + id3;
 
-      int k2 = 1; // Keeps the offset of the permutation relative to k
-      for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
-        if (*it2 != ref_col)
           if (output_type == oJuliaStaticModel)
-            third_derivatives_output << "  @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
-                                     << for_sym.str() << endl;
+            {
+              for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
+              third_derivatives_output << "  @inbounds " << for_sym.str() << " = ";
+              d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
+              third_derivatives_output << endl;
+            }
           else
             {
-              sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
+              sparseHelper(3, third_derivatives_output, k, 0, output_type);
               third_derivatives_output << "=" << eq + 1 << ";" << endl;
 
-              sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
-              third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
+              sparseHelper(3, third_derivatives_output, k, 1, output_type);
+              third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
 
-              sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
-              third_derivatives_output << "=";
               sparseHelper(3, third_derivatives_output, k, 2, output_type);
+              third_derivatives_output << "=";
+              d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
               third_derivatives_output << ";" << endl;
-
-              k2++;
             }
-      k += k2;
+
+          // Compute the column numbers for the 5 other permutations of (id1,id2,id3)
+          // and store them in a set (to avoid duplicates if two indexes are equal)
+          set<int> cols;
+          cols.insert(id1 * hessianColsNbr + id3 * JacobianColsNbr + id2);
+          cols.insert(id2 * hessianColsNbr + id1 * JacobianColsNbr + id3);
+          cols.insert(id2 * hessianColsNbr + id3 * JacobianColsNbr + id1);
+          cols.insert(id3 * hessianColsNbr + id1 * JacobianColsNbr + id2);
+          cols.insert(id3 * hessianColsNbr + id2 * JacobianColsNbr + id1);
+
+          int k2 = 1; // Keeps the offset of the permutation relative to k
+          for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
+            if (*it2 != ref_col)
+              if (output_type == oJuliaStaticModel)
+                third_derivatives_output << "  @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
+                                         << for_sym.str() << endl;
+              else
+                {
+                  sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
+                  third_derivatives_output << "=" << eq + 1 << ";" << endl;
+
+                  sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
+                  third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
+
+                  sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
+                  third_derivatives_output << "=";
+                  sparseHelper(3, third_derivatives_output, k, 2, output_type);
+                  third_derivatives_output << ";" << endl;
+
+                  k2++;
+                }
+          k += k2;
+        }
     }
 
   if (output_type == oMatlabStaticModel)
@@ -1406,69 +1532,94 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
       map<string, string> tmp_paren_vars;
       bool message_printed = false;
       fixNestedParenthesis(model_output, tmp_paren_vars, message_printed);
-      fixNestedParenthesis(model_local_vars_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(model_tt_output, tmp_paren_vars, message_printed);
       fixNestedParenthesis(jacobian_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(jacobian_tt_output, tmp_paren_vars, message_printed);
       fixNestedParenthesis(hessian_output, tmp_paren_vars, message_printed);
+      fixNestedParenthesis(hessian_tt_output, tmp_paren_vars, message_printed);
       fixNestedParenthesis(third_derivatives_output, tmp_paren_vars, message_printed);
-
-      StaticOutput << "residual = zeros( " << equations.size() << ", 1);" << endl << endl
-                   << "%" << endl
-                   << "% Model equations" << endl
-                   << "%" << endl << endl
-                   << model_local_vars_output.str()
-                   << model_output.str()
-                   << "if ~isreal(residual)" << endl
-                   << "  residual = real(residual)+imag(residual).^2;" << endl
-                   << "end" << endl
-                   << "if nargout >= 2," << endl
-                   << "  g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");" << endl << endl
-                   << "  %" << endl
-                   << "  % Jacobian matrix" << endl
-                   << "  %" << endl << endl
-                   << jacobian_output.str()
-                   << "  if ~isreal(g1)" << endl
-                   << "    g1 = real(g1)+2*imag(g1);" << endl
-                   << "  end" << endl
-                   << "if nargout >= 3," << endl
-                   << "  %" << endl
-                   << "  % Hessian matrix" << endl
-                   << "  %" << endl
-                   << endl;
-
+      fixNestedParenthesis(third_derivatives_tt_output, tmp_paren_vars, message_printed);
+
+      string static_name = basename + "_static";
+      ostringstream init_output, end_output;
+      init_output << "residual = zeros(" << equations.size() << ", 1);";
+      end_output << "if ~isreal(residual)" << endl
+                 << "  residual = real(residual)+imag(residual).^2;" << endl
+                 << "end";
+      writeStaticModelHelper(static_name + "_resid", "residual",
+                             static_name + "_resid_tt",
+                             temporary_terms_res_idxs.size(),
+                             "", init_output, end_output,
+                             model_output, model_tt_output);
+
+      init_output.str(string());
+      init_output.clear();
+      end_output.str(string());
+      end_output.clear();
+      init_output << "g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");";
+      end_output << "if ~isreal(g1)" << endl
+                 << "    g1 = real(g1)+2*imag(g1);" << endl
+                 << "end";
+      writeStaticModelHelper(static_name + "_g1", "g1",
+                             static_name + "_g1_tt",
+                             temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size(),
+                             static_name + "_resid_tt",
+                             init_output, end_output,
+                             jacobian_output, jacobian_tt_output);
+      writeWrapperFunctions(static_name, "g1");
+
+      init_output.str(string());
+      init_output.clear();
+      end_output.str(string());
+      end_output.clear();
       if (second_derivatives.size())
-        StaticOutput << "  v2 = zeros(" << NNZDerivatives[1] << ",3);" << endl
-                     << hessian_output.str()
-                     << "  g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << equations.size() << "," << g2ncols << ");" << endl;
+        {
+          init_output << "v2 = zeros(" << NNZDerivatives[1] << ",3);";
+          end_output << "g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << equations.size() << "," << g2ncols << ");";
+        }
       else
-        StaticOutput << "  g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");" << endl;
-
-      // Initialize g3 matrix
-      StaticOutput << "if nargout >= 4," << endl
-                   << "  %" << endl
-                   << "  % Third order derivatives" << endl
-                   << "  %" << endl
-                   << endl;
+        init_output << "g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");";
+      writeStaticModelHelper(static_name + "_g2", "g2",
+                             static_name + "_g2_tt",
+                             temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size()
+                             + temporary_terms_g2_idxs.size(),
+                             static_name + "_g1_tt",
+                             init_output, end_output,
+                             hessian_output, hessian_tt_output);
+      writeWrapperFunctions(static_name, "g2");
+
+      init_output.str(string());
+      init_output.clear();
+      end_output.str(string());
+      end_output.clear();
       int ncols = hessianColsNbr * JacobianColsNbr;
       if (third_derivatives.size())
-        StaticOutput << "  v3 = zeros(" << NNZDerivatives[2] << ",3);" << endl
-                     << third_derivatives_output.str()
-                     << "  g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");" << endl;
-      else // Either 3rd derivatives is all zero, or we didn't compute it
-        StaticOutput << "  g3 = sparse([],[],[]," << nrows << "," << ncols << ");" << endl;
-      StaticOutput << "end" << endl
-                   << "end" << endl
-                   << "end" << endl;
+        {
+          init_output << "v3 = zeros(" << NNZDerivatives[2] << ",3);";
+          end_output << "g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");";
+        }
+      else
+        init_output << "g3 = sparse([],[],[]," << nrows << "," << ncols << ");";
+      writeStaticModelHelper(static_name + "_g3", "g3",
+                             static_name + "_g3_tt",
+                             temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size()
+                             + temporary_terms_g2_idxs.size() + temporary_terms_g3_idxs.size(),
+                             static_name + "_g2_tt",
+                             init_output, end_output,
+                             third_derivatives_output, third_derivatives_tt_output);
+      writeWrapperFunctions(static_name, "g3");
     }
   else if (output_type == oCStaticModel)
     {
+      /*
       StaticOutput << "void Static(double *y, double *x, int nb_row_x, double *params, double *residual, double *g1, double *v2)" << endl
                    << "{" << endl
                    << "  double lhs, rhs;" << endl
                    << endl
-                   << "  /* Residual equations */" << endl
+                   << "  * Residual equations *" << endl
                    << model_local_vars_output.str()
                    << model_output.str()
-                   << "  /* Jacobian  */" << endl
+                   << "  * Jacobian  *" << endl
                    << "  if (g1 == NULL)" << endl
                    << "    return;" << endl
                    << endl
@@ -1476,118 +1627,210 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
                    << endl;
 
       if (second_derivatives.size())
-        StaticOutput << "  /* Hessian for endogenous and exogenous variables */" << endl
+        StaticOutput << "  * Hessian for endogenous and exogenous variables *" << endl
                      << "  if (v2 == NULL)" << endl
                      << "    return;" << endl
                      << endl
                      << hessian_output.str()
                      << endl;
       if (third_derivatives.size())
-        StaticOutput << "  /* Third derivatives for endogenous and exogenous variables */" << endl
+        StaticOutput << "  * Third derivatives for endogenous and exogenous variables *" << endl
                      << "  if (v3 == NULL)" << endl
                      << "    return;" << endl
                      << endl
                      << third_derivatives_output.str()
                      << endl;
+      */
     }
   else
     {
-      ostringstream comments;
-      comments << "## Function Arguments" << endl
-               << endl
-               << "## Input" << endl
-               << " 1 y:         Array{Float64, length(model.endo), 1}             Vector of endogenous variables in declaration order" << endl
-               << " 2 x:         Array{Float64, length(model.exo), 1}              Vector of exogenous variables in declaration order" << endl
-               << " 3 params:    Array{Float64, length(model.param), 1}            Vector of parameter values in declaration order" << endl
-               << endl
-               << "## Output" << endl
-               << " 4 residual:  Array(Float64, model.eq_nbr, 1)                   Vector of residuals of the static model equations" << endl
-               << "                                                                in order of declaration of the equations." << endl
-               << "                                                                Dynare may prepend auxiliary equations, see model.aux_vars" << endl;
-
-      StaticOutput << "function static!(y::Vector{Float64}, x::Vector{Float64}, "
-                   << "params::Vector{Float64}," << endl
-                   << "                 residual::Vector{Float64})" << endl
-                   << "#=" << endl << comments.str() << "=#" << endl
-                   << "  @assert length(y) == " << symbol_table.endo_nbr() << endl
-                   << "  @assert length(x) == " << symbol_table.exo_nbr() << endl
-                   << "  @assert length(params) == " << symbol_table.param_nbr() << endl
-                   << "  @assert length(residual) == " << equations.size() << endl
-                   << "  #" << endl
-                   << "  # Model equations" << endl
-                   << "  #" << endl
-                   << model_local_vars_output.str()
-                   << model_output.str()
-                   << "if ~isreal(residual)" << endl
-                   << "  residual = real(residual)+imag(residual).^2;" << endl
-                   << "end" << endl
-                   << "end" << endl << endl
-                   << "function static!(y::Vector{Float64}, x::Vector{Float64}, "
-                   << "params::Vector{Float64}," << endl
-                   << "                 residual::Vector{Float64}, g1::Matrix{Float64})" << endl;
-
-      comments << " 5 g1:        Array(Float64, model.eq_nbr, length(model.endo))  Jacobian matrix of the static model equations;" << endl
-               << "                                                                columns: variables in declaration order" << endl
-               << "                                                                rows: equations in order of declaration" << endl;
-
-      StaticOutput << "#=" << endl << comments.str() << "=#" << endl
-                   << "  @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr()
-                   << ")" << endl
-                   << "  fill!(g1, 0.0)" << endl
-                   << "  static!(y, x, params, residual)" << endl
-                   << model_local_vars_output.str()
-                   << "  #" << endl
-                   << "  # Jacobian matrix" << endl
-                   << "  #" << endl
-                   << jacobian_output.str()
-                   << "  if ~isreal(g1)" << endl
-                   << "    g1 = real(g1)+2*imag(g1);" << endl
-                   << "  end" << endl
-                   << "end" << endl << endl
-                   << "function static!(y::Vector{Float64}, x::Vector{Float64}, "
-                   << "params::Vector{Float64}," << endl
-                   << "                 residual::Vector{Float64}, g1::Matrix{Float64}, "
-                   << "g2::Matrix{Float64})" << endl;
-
-      comments << " 6 g2:        spzeros(model.eq_nbr, length(model.endo)^2)       Hessian matrix of the static model equations;" << endl
-               << "                                                                columns: variables in declaration order" << endl
-               << "                                                                rows: equations in order of declaration" << endl;
-
-      StaticOutput << "#=" << endl << comments.str() << "=#" << endl
-                   << "  @assert size(g2) == (" << equations.size() << ", " << g2ncols << ")" << endl
-                   << "  fill!(g2, 0.0)" << endl
-                   << "  static!(y, x, params, residual, g1)" << endl;
-      if (second_derivatives.size())
-        StaticOutput << model_local_vars_output.str()
-                     << "  #" << endl
-                     << "  # Hessian matrix" << endl
-                     << "  #" << endl
-                     << hessian_output.str();
+      string filename = basename + "Static.jl";
+      ofstream output;
+      output.open(filename.c_str(), ios::out | ios::binary);
+      if (!output.is_open())
+        {
+          cerr << "Error: Can't open file " << filename << " for writing" << endl;
+          exit(EXIT_FAILURE);
+        }
 
-      // Initialize g3 matrix
+      output << "module " << basename << "Static" << endl
+             << "#" << endl
+             << "# NB: this file was automatically generated by Dynare" << endl
+             << "#     from " << basename << ".mod" << endl
+             << "#" << endl
+             << "using Utils" << endl << endl
+             << "export static!, staticResid!, staticG1!, staticG2!, staticG3!" << endl << endl
+             << "#=" << endl
+             << "# The comments below apply to all functions contained in this module #" << endl
+             << "  NB: The arguments contained on the first line of the function" << endl
+             << "      definition are those that are modified in place" << endl << endl
+             << "## Exported Functions ##" << endl
+             << "  static!      : Wrapper function; computes residuals, Jacobian, Hessian," << endl
+             << "                 and third derivatives depending on the arguments provided" << endl
+             << "  staticResid! : Computes the static model residuals" << endl
+             << "  staticG1!    : Computes the static model Jacobian" << endl
+             << "  staticG2!    : Computes the static model Hessian" << endl
+             << "  staticG3!    : Computes the static model third derivatives" << endl << endl
+             << "## Local Functions ##" << endl
+             << "  staticResidTT! : Computes the static model temporary terms for the residuals" << endl
+             << "  staticG1TT!    : Computes the static model temporary terms for the Jacobian" << endl
+             << "  staticG2TT!    : Computes the static model temporary terms for the Hessian" << endl
+             << "  staticG3TT!    : Computes the static model temporary terms for the third derivatives" << endl << endl
+             << "## Function Arguments ##" << endl
+             << "  T        : Vector{Float64, num_temp_terms, 1}                 Vector of Temporary Terms" << endl
+             << "  y        : Vector{Float64, length(model.endo), 1}             Vector of endogenous variables in declaration order" << endl
+             << "  x        : Vector{Float64, length(model.exo), 1}              Vector of exogenous variables in declaration order" << endl
+             << "  params   : Vector{Float64, length(model.param), 1}            Vector of parameter values in declaration order" << endl
+             << "  residual : Vector{Float64, model.eq_nbr, 1}                   Vector of residuals of the static model equations" << endl
+             << "                                                                in order of declaration of the equations." << endl
+             << "                                                                Dynare may prepend auxiliary equations, see model.aux_vars" << endl
+             << "  g1       : Vector{Float64, model.eq_nbr, length(model.endo)}  Jacobian matrix of the static model equations" << endl
+             << "                                                                columns: variables in declaration order" << endl
+             << "                                                                rows: equations in order of declaration" << endl
+             << "  g2       : spzeros(model.eq_nbr, length(model.endo)^2)        Hessian matrix of the static model equations" << endl
+             << "                                                                columns: variables in declaration order" << endl
+             << "                                                                rows: equations in order of declaration" << endl
+             << "  g3       : spzeros(model.eq_nbr, length(model.endo)^3)        Third derivatives matrix of the static model equations" << endl
+             << "                                                                columns: variables in declaration order" << endl
+             << "                                                                rows: equations in order of declaration" << endl
+             << "=#" << endl << endl;
+
+      // staticResidTT!
+      output << "function staticResidTT!(T::Vector{Float64}," << endl
+             << "                        y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+             << "    @assert length(T) >= " << temporary_terms_res_idxs.size()  << endl
+             << model_tt_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // static!
+      output << "function staticResid!(T::Vector{Float64}, residual::Vector{Float64}," << endl
+             << "                      y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T_flag::Bool)" << endl
+             << "    @assert length(y) == " << symbol_table.endo_nbr() << endl
+             << "    @assert length(x) == " << symbol_table.exo_nbr() << endl
+             << "    @assert length(params) == " << symbol_table.param_nbr() << endl
+             << "    @assert length(residual) == " << equations.size() << endl
+             << "    if T_flag" << endl
+             << "        staticResidTT!(T, y, x, params)" << endl
+             << "    end" << endl
+             << model_output.str()
+             << "    if ~isreal(residual)" << endl
+             << "        residual = real(residual)+imag(residual).^2;" << endl
+             << "    end" << endl
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // staticG1TT!
+      output << "function staticG1TT!(T::Vector{Float64}," << endl
+             << "                     y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+             << "    staticResidTT!(T, y, x, params)" << endl
+             << jacobian_tt_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // staticG1!
+      output << "function staticG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl
+             << "                   y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T_flag::Bool)" << endl
+             << "    @assert length(T) >= "
+             << temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size() << endl
+             << "    @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr() << ")" << endl
+             << "    @assert length(y) == " << symbol_table.endo_nbr() << endl
+             << "    @assert length(x) == " << symbol_table.exo_nbr() << endl
+             << "    @assert length(params) == " << symbol_table.param_nbr() << endl
+             << "    if T_flag" << endl
+             << "        staticG1TT!(T, y, x, params)" << endl
+             << "    end" << endl
+             << "    fill!(g1, 0.0)" << endl
+             << jacobian_output.str()
+             << "    if ~isreal(g1)" << endl
+             << "        g1 = real(g1)+2*imag(g1);" << endl
+             << "    end" << endl
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // staticG2TT!
+      output << "function staticG2TT!(T::Vector{Float64}," << endl
+             << "                     y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+             << "    staticG1TT!(T, y, x, params)" << endl
+             << hessian_tt_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // staticG2!
+      output << "function staticG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl
+             << "                   y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T_flag::Bool)" << endl
+             << "    @assert length(T) >= "
+             << temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size() + temporary_terms_g2_idxs.size() << endl
+             << "    @assert size(g2) == (" << equations.size() << ", " << g2ncols << ")" << endl
+             << "    @assert length(y) == " << symbol_table.endo_nbr() << endl
+             << "    @assert length(x) == " << symbol_table.exo_nbr() << endl
+             << "    @assert length(params) == " << symbol_table.param_nbr() << endl
+             << "    if T_flag" << endl
+             << "        staticG2TT!(T, y, x, params)" << endl
+             << "    end" << endl
+             << "    fill!(g2, 0.0)" << endl
+             << hessian_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // staticG3TT!
+      output << "function staticG3TT!(T::Vector{Float64}," << endl
+             << "                     y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+             << "    staticG2TT!(T, y, x, params)" << endl
+             << third_derivatives_tt_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // staticG3!
       int ncols = hessianColsNbr * JacobianColsNbr;
-      StaticOutput << "end" << endl << endl
-                   << "function static!(y::Vector{Float64}, x::Vector{Float64}, "
-                   << "params::Vector{Float64}," << endl
-                   << "                 residual::Vector{Float64}, g1::Matrix{Float64}, "
-                   << "g2::Matrix{Float64}," << endl
-                   << "                 g3::Matrix{Float64})" << endl;
-
-      comments << " 7 g3:        spzeros(model.eq_nbr, length(model.endo)^3)       Third derivatives matrix of the static model equations;" << endl
-               << "                                                                columns: variables in declaration order" << endl
-               << "                                                                rows: equations in order of declaration" << endl;
-
-      StaticOutput << "#=" << endl << comments.str() << "=#" << endl
-                   << "  @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
-                   << "  fill!(g3, 0.0)" << endl
-                   << "  static!(y, x, params, residual, g1, g2)" << endl;
-      if (third_derivatives.size())
-        StaticOutput << model_local_vars_output.str()
-                     << "  #" << endl
-                     << "  # Third order derivatives" << endl
-                     << "  #" << endl
-                     << third_derivatives_output.str();
-      StaticOutput << "end" << endl;
+      output << "function staticG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl
+             << "                   y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T_flag::Bool)" << endl
+             << "    @assert length(T) >= "
+             << temporary_terms_res_idxs.size() + temporary_terms_g1_idxs.size() + temporary_terms_g2_idxs.size() + temporary_terms_g3_idxs.size() << endl
+             << "    @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
+             << "    @assert length(y) == " << symbol_table.endo_nbr() << endl
+             << "    @assert length(x) == " << symbol_table.exo_nbr() << endl
+             << "    @assert length(params) == " << symbol_table.param_nbr() << endl
+             << "    if T_flag" << endl
+             << "        staticG3TT!(T, y, x, params)" << endl
+             << "    end" << endl
+             << "    fill!(g3, 0.0)" << endl
+             << third_derivatives_output.str()
+             << "    return nothing" << endl
+             << "end" << endl << endl;
+
+      // static!
+      output  << "function static!(T::Vector{Float64}, residual::Vector{Float64}," << endl
+             << "                  y::Vector{Float64}, x::Matrix{Float64}, params::Vector{Float64})" << endl
+             << "    staticResid!(T, residual, y, x, params, true)" << endl
+             << "    return nothing" << endl
+             << "end" << endl
+             << endl
+             << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}," << endl
+             << "                 y::Vector{Float64}, x::Matrix{Float64}, params::Vector{Float64})" << endl
+             << "    staticG1!(T, g1, y, x, params, true)" << endl
+             << "    staticResid!(T, residual, y, x, params, false)" << endl
+             << "    return nothing" << endl
+             << "end" << endl
+             << endl
+             << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl
+             << "                 y::Vector{Float64}, x::Matrix{Float64}, params::Vector{Float64})" << endl
+             << "    staticG2!(T, g2, y, x, params, true)" << endl
+             << "    staticG1!(T, g1, y, x, params, false)" << endl
+             << "    staticResid!(T, residual, y, x, params, false)" << endl
+             << "    return nothing" << endl
+             << "end" << endl
+             << endl
+             << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}, g3::Matrix{Float64}," << endl
+             << "                 y::Vector{Float64}, x::Matrix{Float64}, params::Vector{Float64})" << endl
+             << "    staticG3!(T, g3, y, x, params, true)" << endl
+             << "    staticG2!(T, g2, y, x, params, false)" << endl
+             << "    staticG1!(T, g1, y, x, params, false)" << endl
+             << "    staticResid!(T, residual, y, x, params, false)" << endl
+             << "    return nothing" << endl
+             << "end" << endl
+             << "end" << endl;
+      output.close();
     }
 }
 
@@ -1711,24 +1954,7 @@ StaticModel::writeStaticCFile(const string &func_name) const
 void
 StaticModel::writeStaticJuliaFile(const string &basename) const
 {
-  string filename = basename + "Static.jl";
-  ofstream output;
-  output.open(filename.c_str(), ios::out | ios::binary);
-  if (!output.is_open())
-    {
-      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  output << "module " << basename << "Static" << endl
-         << "#" << endl
-         << "# NB: this file was automatically generated by Dynare" << endl
-         << "#     from " << basename << ".mod" << endl
-         << "#" << endl
-         << "using Utils" << endl << endl
-         << "export static!" << endl << endl;
-  writeStaticModel(output, false, true);
-  output << "end" << endl;
+  writeStaticModel(basename, false, true);
 }
 
 void
@@ -2171,10 +2397,11 @@ StaticModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType
 {
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temporary_terms;
+  temporary_terms_idxs_t temporary_terms_idxs;
   for (int i = 0; i < (int) aux_equations.size(); i++)
     if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
       dynamic_cast<ExprNode *>(aux_equations[i])->writeExternalFunctionOutput(output, oMatlabStaticModel,
-                                                                              temporary_terms, tef_terms);
+                                                                              temporary_terms, temporary_terms_idxs, tef_terms);
   for (int i = 0; i < (int) aux_equations.size(); i++)
     {
       dynamic_cast<ExprNode *>(aux_equations[i]->substituteStaticAuxiliaryDefinition())->writeOutput(output, output_type, temporary_terms, tef_terms);
@@ -2187,10 +2414,11 @@ StaticModel::writeLatexAuxVarRecursiveDefinitions(ostream &output) const
 {
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temporary_terms;
+    temporary_terms_idxs_t temporary_terms_idxs;
   for (int i = 0; i < (int) aux_equations.size(); i++)
     if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
       dynamic_cast<ExprNode *>(aux_equations[i])->writeExternalFunctionOutput(output, oLatexStaticModel,
-                                                                              temporary_terms, tef_terms);
+                                                                              temporary_terms, temporary_terms_idxs, tef_terms);
   for (int i = 0; i < (int) aux_equations.size(); i++)
     {
       output << "\\begin{dmath}" << endl;
@@ -2490,7 +2718,6 @@ StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) co
   ostringstream third_derivatives_output;  // Used for storing third order derivatives equations
 
   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;
 
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index ad55f89d90c2e2b188d1ef4d2fd7dcc9b510990a..54394d0e353dad77392d95377d0d4082b6e43318 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -51,7 +51,7 @@ private:
   void writeStaticJuliaFile(const string &basename) const;
 
   //! Writes the static model equations and its derivatives
-  void writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const;
+  void writeStaticModel(const string &basename, ostream &StaticOutput, bool use_dll, bool julia) const;
 
   //! Writes the static function calling the block to solve (Matlab version)
   void writeStaticBlockMFSFile(const string &basename) const;
@@ -148,6 +148,15 @@ protected:
   //!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
   vector<pair<int, int> > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
 
+  //! Helper functions for writeStaticModel
+  void writeStaticModelHelper(const string &name, const string &retvalname,
+                              const string &name_tt, size_t ttlen,
+                              const string &previous_tt_name,
+                              const ostringstream &init_s, const ostringstream &end_s,
+                              const ostringstream &s, const ostringstream &s_tt) const;
+  void writeWrapperFunctions(const string &basename, const string &ending) const;
+  void writeStaticModel(ostream &DynamicOutput, bool use_dll, bool julia) const;
+  void writeStaticModel(const string &dynamic_basename, bool use_dll, bool julia) const;
 public:
   StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants, ExternalFunctionsTable &external_functions_table_arg);