diff --git a/DynamicModel.cc b/DynamicModel.cc
index 00cb992544f1a5f9b6c744ff7e5c48911dbc65b4..7dfc7d7ecb305dcacb336b65e6bfbda0ae5e7095 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -4284,281 +4284,4 @@ DynamicModel::dynamicOnlyEquationsNbr() const
   return eqs.size();
 }
 
-void
-DynamicModel::writeFirstDerivativesC(const string &basename, bool cuda) const
-{
-  string filename = basename + "_first_derivatives.c";
-  ofstream mDynamicModelFile, mDynamicMexFile;
-
-  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 << "/*" << endl
-                    << " * " << filename << " : Computes first order derivatives of the model for Dynare" << endl
-                    << " *" << endl
-                    << " * Warning : this file is generated automatically by Dynare" << endl
-                    << " *           from model " << basename << "(.mod)" << endl
-                    << " */" << endl
-                    << "#include <math.h>" << endl;
-
-  mDynamicModelFile << "#include <stdlib.h>" << endl;
-
-  mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
-                    << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
-
-  // Write function definition if oPowerDeriv is used
-  writePowerDerivCHeader(mDynamicModelFile);
-
-  mDynamicModelFile << "void FirstDerivatives(const 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;
-
-  // this is always empty here, but needed by d1->writeOutput
-  deriv_node_temp_terms_t tef_terms;
-
-  // Writing Jacobian
-  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(mDynamicModelFile, eq, getDynJacobianCol(var), oCDynamicModel);
-      mDynamicModelFile << "=";
-      // oCstaticModel makes reference to the static variables
-      d1->writeOutput(mDynamicModelFile, oCStaticModel, temporary_terms, tef_terms);
-      mDynamicModelFile << ";" << endl;
-    }
-  
-  mDynamicModelFile << "}" << endl;
-
-  writePowerDeriv(mDynamicModelFile, true);
-  mDynamicModelFile.close();
-
-}
-
-// using compressed sparse row format (CSR)
-void
-DynamicModel::writeSecondDerivativesC_csr(const string &basename, bool cuda) const
-{
-
-  string filename = basename + "_second_derivatives.c";
-  ofstream mDynamicModelFile, mDynamicMexFile;
-
-  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 << "/*" << endl
-                    << " * " << filename << " : Computes second order derivatives of the model for Dynare" << endl
-                    << " *" << endl
-                    << " * Warning : this file is generated automatically by Dynare" << endl
-                    << " *           from model " << basename << "(.mod)" << endl
-                    << " */" << endl
-                    << "#include <math.h>" << endl;
-
-  mDynamicModelFile << "#include <stdlib.h>" << endl;
-
-  mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
-                    << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
-
-  // write function definition if oPowerDeriv is used
-  writePowerDerivCHeader(mDynamicModelFile);
-
-  mDynamicModelFile << "void SecondDerivatives(const double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, int *row_ptr, int *col_ptr, double *value)" << endl
-                    << "{" << endl;
-
-  // this is always empty here, but needed by d1->writeOutput
-  deriv_node_temp_terms_t tef_terms;
-
-  // Indexing derivatives in column order
-  vector<derivative> D;
-  int hessianColsNbr = dynJacobianColsNbr*dynJacobianColsNbr;
-  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;
-
-      int id1 = getDynJacobianCol(var1);
-      int id2 = getDynJacobianCol(var2);
-
-      int col_nb = id1 * dynJacobianColsNbr + id2;
-
-      derivative deriv(col_nb + eq*hessianColsNbr,col_nb,eq,it->second);
-      D.push_back(deriv);
-      if (id1 != id2)
-	{
-	  col_nb = id2 * dynJacobianColsNbr + id1;
-	  derivative deriv(col_nb + eq*hessianColsNbr,col_nb,eq,it->second);
-	  D.push_back(deriv);
-	}
-    }
-  sort(D.begin(), D.end(), derivative_less_than() );
-
-  // Writing Hessian
-  vector<int> row_ptr(equations.size());
-  fill(row_ptr.begin(),row_ptr.end(),0.0);
-  int k = 0;
-  for(vector<derivative>::const_iterator it = D.begin(); it != D.end(); ++it)
-    {
-      row_ptr[it->row_nbr]++;
-      mDynamicModelFile << "col_ptr[" << k << "] "
-			<< "=" << it->col_nbr << ";" << endl;
-      mDynamicModelFile << "value[" << k << "] = ";
-      // oCstaticModel makes reference to the static variables
-      it->value->writeOutput(mDynamicModelFile, oCStaticModel, temporary_terms, tef_terms);
-      mDynamicModelFile << ";" << endl;
-      k++;
-    }
-  
-  // row_ptr must point to the relative address of the first element of the row
-  int cumsum = 0;
-  mDynamicModelFile << "row_ptr = [ 0";
-  for (vector<int>::iterator it=row_ptr.begin(); it != row_ptr.end(); ++it)
-    {
-      cumsum += *it;
-      mDynamicModelFile << ", " << cumsum;
-    }
-  mDynamicModelFile << "];" << endl;   
-
-  mDynamicModelFile << "}" << endl;
-
-  writePowerDeriv(mDynamicModelFile, true);
-  mDynamicModelFile.close();
-
-}
-
-void
-DynamicModel::writeThirdDerivativesC_csr(const string &basename, bool cuda) const
-{
-  string filename = basename + "_third_derivatives.c";
-  ofstream mDynamicModelFile, mDynamicMexFile;
-
-  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 << "/*" << endl
-                    << " * " << filename << " : Computes third order derivatives of the model for Dynare" << endl
-                    << " *" << endl
-                    << " * Warning : this file is generated automatically by Dynare" << endl
-                    << " *           from model " << basename << "(.mod)" << endl
-                    << " */" << endl
-                    << "#include <math.h>" << endl;
-
-  mDynamicModelFile << "#include <stdlib.h>" << endl;
-
-  mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
-                    << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
-
-  // Write function definition if oPowerDeriv is used
-  writePowerDerivCHeader(mDynamicModelFile);
-
-  mDynamicModelFile << "void ThirdDerivatives(const 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;
-
-  // this is always empty here, but needed by d1->writeOutput
-  deriv_node_temp_terms_t tef_terms;
-
-  vector<derivative> D;
-  int hessianColsNbr = dynJacobianColsNbr*dynJacobianColsNbr;
-  int thirdDerivativesColsNbr = hessianColsNbr*dynJacobianColsNbr;
-  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;
-
-      int id1 = getDynJacobianCol(var1);
-      int id2 = getDynJacobianCol(var2);
-      int id3 = getDynJacobianCol(var3);
-
-      // Reference column number for the g3 matrix (with symmetrical derivatives)
-      vector<long unsigned int>  cols;
-      long unsigned int col_nb = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
-      int thirdDColsNbr = hessianColsNbr*dynJacobianColsNbr;
-      derivative deriv(col_nb + eq*thirdDColsNbr,col_nb,eq,it->second);
-      D.push_back(deriv);
-      cols.push_back(col_nb);
-      col_nb = id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2;
-      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
-	{
-	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
-	  D.push_back(deriv);
-	  cols.push_back(col_nb);
-	}
-      col_nb = id2 * hessianColsNbr + id1 * dynJacobianColsNbr + id3;
-      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
-	{
-	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
-	  D.push_back(deriv);
-	  cols.push_back(col_nb);
-	}
-      col_nb = id2 * hessianColsNbr + id3 * dynJacobianColsNbr + id1;
-      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
-	{
-	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
-	  D.push_back(deriv);
-	  cols.push_back(col_nb);
-	}
-      col_nb = id3 * hessianColsNbr + id1 * dynJacobianColsNbr + id2;
-      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
-	{
-	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
-	  D.push_back(deriv);
-	  cols.push_back(col_nb);
-	}
-      col_nb = id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1;
-      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
-	{
-	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
-	  D.push_back(deriv);
-	}
-    }
-
-  sort(D.begin(), D.end(), derivative_less_than() );
-
-  vector<int> row_ptr(equations.size());
-  fill(row_ptr.begin(),row_ptr.end(),0.0);
-  int k = 0;
-  for(vector<derivative>::const_iterator it = D.begin(); it != D.end(); ++it)
-    {
-      row_ptr[it->row_nbr]++;
-      mDynamicModelFile << "col_ptr[" << k << "] "
-			<< "=" << it->col_nbr << ";" << endl;
-      mDynamicModelFile << "value[" << k << "] = ";
-      // oCstaticModel makes reference to the static variables
-      it->value->writeOutput(mDynamicModelFile, oCStaticModel, temporary_terms, tef_terms);
-      mDynamicModelFile << ";" << endl;
-      k++;
-    }
-
-  // row_ptr must point to the relative address of the first element of the row
-  int cumsum = 0;
-  mDynamicModelFile << "row_ptr = [ 0";
-  for (vector<int>::iterator it=row_ptr.begin(); it != row_ptr.end(); ++it)
-    {
-      cumsum += *it;
-      mDynamicModelFile << ", " << cumsum;
-    }
-  mDynamicModelFile << "];" << endl;   
-
-  mDynamicModelFile << "}" << endl;
-
-  writePowerDeriv(mDynamicModelFile, true);
-  mDynamicModelFile.close();
-
-}
-
 
diff --git a/DynamicModel.hh b/DynamicModel.hh
index e1ee490bb130d91afa51119ec5e9d03bdf7e9c4b..c752fe0bcf4855f687703e37c2ee21acefaa7e61 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -213,10 +213,6 @@ public:
                      const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode);
   //! Writes model initialization and lead/lag incidence matrix to output
   void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const;
-  //! Writes model initialization and lead/lag incidence matrix to C output
-  void writeCOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const;
-  //! Writes model initialization and lead/lag incidence matrix to Cpp output
-  void writeCCOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const;
 
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const string &bin_basename,
@@ -225,12 +221,6 @@ public:
   void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order) const;
   //! Writes file containing parameters derivatives
   void writeParamsDerivativesFile(const string &basename) const;
-  //! Writes CC file containing first order derivatives of model evaluated at steady state
-  void writeFirstDerivativesC(const string &basename, bool cuda) const;
-  //! Writes CC file containing second order derivatives of model evaluated at steady state (compressed sparse column)
-  void writeSecondDerivativesC_csr(const string &basename, bool cuda) const;
-  //! Writes CC file containing third order derivatives of model evaluated at steady state (compressed sparse column)
-  void writeThirdDerivativesC_csr(const string &basename, bool cuda) const;
   //! Converts to static model (only the equations)
   /*! It assumes that the static model given in argument has just been allocated */
   void toStatic(StaticModel &static_model) const;
@@ -473,6 +463,18 @@ public:
       return (-1);
   };
   bool isModelLocalVariableUsed() const;
+
+  // in ExternalFiles.cc
+  //! Writes model initialization and lead/lag incidence matrix to C output
+  void writeCOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const;
+  //! Writes model initialization and lead/lag incidence matrix to Cpp output
+  void writeCCOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const;
+  //! Writes C file containing first order derivatives of model evaluated at steady state
+  void writeFirstDerivativesC(const string &basename, bool cuda) const;
+  //! Writes C file containing second order derivatives of model evaluated at steady state (compressed sparse column)
+  void writeSecondDerivativesC_csr(const string &basename, bool cuda) const;
+  //! Writes C file containing third order derivatives of model evaluated at steady state (compressed sparse column)
+  void writeThirdDerivativesC_csr(const string &basename, bool cuda) const;
 };
 
 //! Classes to re-order derivatives for various sparse storage formats 
diff --git a/ExternalFiles.cc b/ExternalFiles.cc
index 76dc3d30bf6a28b059bce0e30ee5e1f5d3b58508..9ac2ba72056869333e9dfd4b2c0bc5ed648b5ca8 100644
--- a/ExternalFiles.cc
+++ b/ExternalFiles.cc
@@ -16,6 +16,7 @@
  * You should have received a copy of the GNU General Public License
  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
  */
+#include <cassert>
 #include "ModFile.hh"
 #include "DynamicModel.hh"
 #include "StaticModel.hh"
@@ -154,14 +155,13 @@ DynamicModel::writeCOutput(ostream &output, const string &basename, bool block_d
   vector<int> zeta_back, zeta_mixed, zeta_fwrd, zeta_static;
   for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++)
     {
-      int varID;
       // Loop on periods
       for (int lag = 0; lag <= 2; lag++)
 	{
 	  lag_presence[lag] = 1;
           try
             {
-              varID = getDerivID(symbol_table.getID(eEndogenous, endoID), lag-1);
+              getDerivID(symbol_table.getID(eEndogenous, endoID), lag-1);
             }
           catch (UnknownDerivIDException &e)
             {
@@ -239,7 +239,335 @@ DynamicModel::writeCOutput(ostream &output, const string &basename, bool block_d
     }
 }
 
+void
+DynamicModel::writeFirstDerivativesC(const string &basename, bool cuda) const
+{
+  string filename = basename + "_first_derivatives.c";
+  ofstream mDynamicModelFile, mDynamicMexFile;
+
+  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 << "/*" << endl
+                    << " * " << filename << " : Computes first order derivatives of the model for Dynare" << endl
+                    << " *" << endl
+                    << " * Warning : this file is generated automatically by Dynare" << endl
+                    << " *           from model " << basename << "(.mod)" << endl
+                    << " */" << endl
+                    << "#include <math.h>" << endl;
+
+  mDynamicModelFile << "#include <stdlib.h>" << endl;
+
+  mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
+                    << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
+
+  // Write function definition if oPowerDeriv is used
+  writePowerDerivCHeader(mDynamicModelFile);
+
+  mDynamicModelFile << "void FirstDerivatives(const 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;
+
+  // this is always empty here, but needed by d1->writeOutput
+  deriv_node_temp_terms_t tef_terms;
+
+  // Writing Jacobian
+  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(mDynamicModelFile, eq, getDynJacobianCol(var), oCDynamicModel);
+      mDynamicModelFile << "=";
+      // oCstaticModel makes reference to the static variables
+      d1->writeOutput(mDynamicModelFile, oCStaticModel, temporary_terms, tef_terms);
+      mDynamicModelFile << ";" << endl;
+    }
+  
+  mDynamicModelFile << "}" << endl;
+
+  writePowerDeriv(mDynamicModelFile, true);
+  mDynamicModelFile.close();
+
+}
+
+// using compressed sparse row format (CSR)
+void
+DynamicModel::writeSecondDerivativesC_csr(const string &basename, bool cuda) const
+{
+
+  string filename = basename + "_second_derivatives.c";
+  ofstream mDynamicModelFile, mDynamicMexFile;
+
+  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 << "/*" << endl
+                    << " * " << filename << " : Computes second order derivatives of the model for Dynare" << endl
+                    << " *" << endl
+                    << " * Warning : this file is generated automatically by Dynare" << endl
+                    << " *           from model " << basename << "(.mod)" << endl
+                    << " */" << endl
+                    << "#include <math.h>" << endl;
+
+  mDynamicModelFile << "#include <stdlib.h>" << endl;
+
+  mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
+                    << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
+
+  // write function definition if oPowerDeriv is used
+  writePowerDerivCHeader(mDynamicModelFile);
+
+  mDynamicModelFile << "void SecondDerivatives(const double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, int *row_ptr, int *col_ptr, double *value)" << endl
+                    << "{" << endl;
+
+  // this is always empty here, but needed by d1->writeOutput
+  deriv_node_temp_terms_t tef_terms;
+
+  // Indexing derivatives in column order
+  vector<derivative> D;
+  int hessianColsNbr = dynJacobianColsNbr*dynJacobianColsNbr;
+  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;
+
+      int id1 = getDynJacobianCol(var1);
+      int id2 = getDynJacobianCol(var2);
+
+      int col_nb = id1 * dynJacobianColsNbr + id2;
+
+      derivative deriv(col_nb + eq*hessianColsNbr,col_nb,eq,it->second);
+      D.push_back(deriv);
+      if (id1 != id2)
+	{
+	  col_nb = id2 * dynJacobianColsNbr + id1;
+	  derivative deriv(col_nb + eq*hessianColsNbr,col_nb,eq,it->second);
+	  D.push_back(deriv);
+	}
+    }
+  sort(D.begin(), D.end(), derivative_less_than() );
+
+  // Writing Hessian
+  vector<int> row_ptr(equations.size());
+  fill(row_ptr.begin(),row_ptr.end(),0.0);
+  int k = 0;
+  for(vector<derivative>::const_iterator it = D.begin(); it != D.end(); ++it)
+    {
+      row_ptr[it->row_nbr]++;
+      mDynamicModelFile << "col_ptr[" << k << "] "
+			<< "=" << it->col_nbr << ";" << endl;
+      mDynamicModelFile << "value[" << k << "] = ";
+      // oCstaticModel makes reference to the static variables
+      it->value->writeOutput(mDynamicModelFile, oCStaticModel, temporary_terms, tef_terms);
+      mDynamicModelFile << ";" << endl;
+      k++;
+    }
+  
+  // row_ptr must point to the relative address of the first element of the row
+  int cumsum = 0;
+  mDynamicModelFile << "row_ptr = [ 0";
+  for (vector<int>::iterator it=row_ptr.begin(); it != row_ptr.end(); ++it)
+    {
+      cumsum += *it;
+      mDynamicModelFile << ", " << cumsum;
+    }
+  mDynamicModelFile << "];" << endl;   
+
+  mDynamicModelFile << "}" << endl;
+
+  writePowerDeriv(mDynamicModelFile, true);
+  mDynamicModelFile.close();
+
+}
+
+void
+DynamicModel::writeThirdDerivativesC_csr(const string &basename, bool cuda) const
+{
+  string filename = basename + "_third_derivatives.c";
+  ofstream mDynamicModelFile, mDynamicMexFile;
+
+  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 << "/*" << endl
+                    << " * " << filename << " : Computes third order derivatives of the model for Dynare" << endl
+                    << " *" << endl
+                    << " * Warning : this file is generated automatically by Dynare" << endl
+                    << " *           from model " << basename << "(.mod)" << endl
+                    << " */" << endl
+                    << "#include <math.h>" << endl;
+
+  mDynamicModelFile << "#include <stdlib.h>" << endl;
+
+  mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
+                    << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
+
+  // Write function definition if oPowerDeriv is used
+  writePowerDerivCHeader(mDynamicModelFile);
+
+  mDynamicModelFile << "void ThirdDerivatives(const 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;
+
+  // this is always empty here, but needed by d1->writeOutput
+  deriv_node_temp_terms_t tef_terms;
+
+  vector<derivative> D;
+  int hessianColsNbr = dynJacobianColsNbr*dynJacobianColsNbr;
+  int thirdDerivativesColsNbr = hessianColsNbr*dynJacobianColsNbr;
+  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;
+
+      int id1 = getDynJacobianCol(var1);
+      int id2 = getDynJacobianCol(var2);
+      int id3 = getDynJacobianCol(var3);
+
+      // Reference column number for the g3 matrix (with symmetrical derivatives)
+      vector<long unsigned int>  cols;
+      long unsigned int col_nb = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
+      int thirdDColsNbr = hessianColsNbr*dynJacobianColsNbr;
+      derivative deriv(col_nb + eq*thirdDColsNbr,col_nb,eq,it->second);
+      D.push_back(deriv);
+      cols.push_back(col_nb);
+      col_nb = id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2;
+      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
+	{
+	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
+	  D.push_back(deriv);
+	  cols.push_back(col_nb);
+	}
+      col_nb = id2 * hessianColsNbr + id1 * dynJacobianColsNbr + id3;
+      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
+	{
+	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
+	  D.push_back(deriv);
+	  cols.push_back(col_nb);
+	}
+      col_nb = id2 * hessianColsNbr + id3 * dynJacobianColsNbr + id1;
+      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
+	{
+	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
+	  D.push_back(deriv);
+	  cols.push_back(col_nb);
+	}
+      col_nb = id3 * hessianColsNbr + id1 * dynJacobianColsNbr + id2;
+      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
+	{
+	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
+	  D.push_back(deriv);
+	  cols.push_back(col_nb);
+	}
+      col_nb = id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1;
+      if (find(cols.begin(),cols.end(),col_nb) == cols.end())
+	{
+	  derivative deriv(col_nb + eq*thirdDerivativesColsNbr,col_nb,eq,it->second);
+	  D.push_back(deriv);
+	}
+    }
+
+  sort(D.begin(), D.end(), derivative_less_than() );
+
+  vector<int> row_ptr(equations.size());
+  fill(row_ptr.begin(),row_ptr.end(),0.0);
+  int k = 0;
+  for(vector<derivative>::const_iterator it = D.begin(); it != D.end(); ++it)
+    {
+      row_ptr[it->row_nbr]++;
+      mDynamicModelFile << "col_ptr[" << k << "] "
+			<< "=" << it->col_nbr << ";" << endl;
+      mDynamicModelFile << "value[" << k << "] = ";
+      // oCstaticModel makes reference to the static variables
+      it->value->writeOutput(mDynamicModelFile, oCStaticModel, temporary_terms, tef_terms);
+      mDynamicModelFile << ";" << endl;
+      k++;
+    }
+
+  // row_ptr must point to the relative address of the first element of the row
+  int cumsum = 0;
+  mDynamicModelFile << "row_ptr = [ 0";
+  for (vector<int>::iterator it=row_ptr.begin(); it != row_ptr.end(); ++it)
+    {
+      cumsum += *it;
+      mDynamicModelFile << ", " << cumsum;
+    }
+  mDynamicModelFile << "];" << endl;   
+
+  mDynamicModelFile << "}" << endl;
+
+  writePowerDeriv(mDynamicModelFile, true);
+  mDynamicModelFile.close();
+
+}
+
+void
+SteadyStateModel::writeSteadyStateFileC(const string &basename, bool ramsey_model) const
+{
+  string filename = basename + "_steadystate.c";
+
+  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 << "#include <math.h>" << endl;
+
+  output << "void steadystate("
+	 << "const double *exo_, const double *params, double *ys_, int *info)" << endl
+         << "// Steady state file generated by Dynare preprocessor" << endl
+	 << "{" << endl
+         << "    *info = 0;" << endl;
+
+  if (def_table.size() == 0)
+    {
+      output << "    return;" << endl
+	     << "}" << endl;
+      return;
+    }
+
+  for (size_t i = 0; i < def_table.size(); i++)
+    {
+      const vector<int> &symb_ids = def_table[i].first;
+      output << "    ";
+      if (symb_ids.size() > 1)
+	std::cout << "Error: in C, multiple returns are not permitted in steady_state_model" << std::endl;
+      variable_node_map_t::const_iterator it = variable_node_map.find(make_pair(symb_ids[0], 0));
+      assert(it != variable_node_map.end());
+      if (it->second->get_type() == eModFileLocalVariable)
+	output << "double ";
+      dynamic_cast<ExprNode *>(it->second)->writeOutput(output, oCSteadyStateFile);
+      output << "=";
+      def_table[i].second->writeOutput(output, oCSteadyStateFile);
+      output << ";" << endl;
+    }
+  output << "    // Auxiliary equations" << endl;
+  static_model.writeAuxVarInitval(output, oCSteadyStateFile);
+  output << "}" << endl;
+}
+
+
+//
 // C++ interface
+//
 void
 ModFile::writeExternalFilesCC(const string &basename, FileOutputType output) const
 {
@@ -353,14 +681,13 @@ DynamicModel::writeCCOutput(ostream &output, const string &basename, bool block_
   // Loop on endogenous variables
   for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++)
     {
-      int varID;
       // Loop on periods
       for (int lag = 0; lag <= 2; lag++)
 	{
 	  lag_presence[lag] = 1;
           try
             {
-              varID = getDerivID(symbol_table.getID(eEndogenous, endoID), lag-1);
+              getDerivID(symbol_table.getID(eEndogenous, endoID), lag-1);
             }
           catch (UnknownDerivIDException &e)
             {
diff --git a/SteadyStateModel.cc b/SteadyStateModel.cc
index 1c811c523f5aa0170d05439f1e77546694413b1f..d19a7f38c8a014c292ea35591c86b75ca260b7d3 100644
--- a/SteadyStateModel.cc
+++ b/SteadyStateModel.cc
@@ -152,51 +152,3 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
          << "end" << endl;
 }
 
-void
-SteadyStateModel::writeSteadyStateFileC(const string &basename, bool ramsey_model) const
-{
-  string filename = basename + "_steadystate.c";
-
-  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 << "#include <math.h>" << endl;
-
-  output << "void steadystate("
-	 << "const double *exo_, const double *params, double *ys_, int *info)" << endl
-         << "// Steady state file generated by Dynare preprocessor" << endl
-	 << "{" << endl
-         << "    *info = 0;" << endl;
-
-  if (def_table.size() == 0)
-    {
-      output << "    return;" << endl
-	     << "}" << endl;
-      return;
-    }
-
-  for (size_t i = 0; i < def_table.size(); i++)
-    {
-      const vector<int> &symb_ids = def_table[i].first;
-      output << "    ";
-      if (symb_ids.size() > 1)
-	std::cout << "Error: in C, multiple returns are not permitted in steady_state_model" << std::endl;
-      variable_node_map_t::const_iterator it = variable_node_map.find(make_pair(symb_ids[0], 0));
-      assert(it != variable_node_map.end());
-      if (it->second->get_type() == eModFileLocalVariable)
-	output << "double ";
-      dynamic_cast<ExprNode *>(it->second)->writeOutput(output, oCSteadyStateFile);
-      output << "=";
-      def_table[i].second->writeOutput(output, oCSteadyStateFile);
-      output << ";" << endl;
-    }
-  output << "    // Auxiliary equations" << endl;
-  static_model.writeAuxVarInitval(output, oCSteadyStateFile);
-  output << "}" << endl;
-}
-
diff --git a/SteadyStateModel.hh b/SteadyStateModel.hh
index 9a28c430e069867c82285d9e8fff324dc88185e5..64b78573b88ae2ae39affa803cb372e3dc1c1b64 100644
--- a/SteadyStateModel.hh
+++ b/SteadyStateModel.hh
@@ -49,6 +49,7 @@ public:
     \param[in] ramsey_model Is there a Ramsey model in the MOD file? If yes, then use the "ys" in argument of the steady state file as initial values
   */
   void writeSteadyStateFile(const string &basename, bool ramsey_model) const;
+  // in ExternalFiles.cc
   void writeSteadyStateFileC(const string &basename, bool ramsey_model) const;
 };