diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 8614dda5d119cc98962ba2d8d7c6b617b7c8b674..cf12f32965fa8334d057a44bad78ca9d6d073bdb 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -5750,587 +5750,6 @@ DynamicModel::isChecksumMatching(const string &basename) const
   return true;
 }
 
-void
-DynamicModel::writeCOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present) const
-{
-  int lag_presence[3];
-  // Loop on endogenous variables
-  vector<int> zeta_back, zeta_mixed, zeta_fwrd, zeta_static;
-  for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++)
-    {
-      // Loop on periods
-      for (int lag = 0; lag <= 2; lag++)
-        {
-          lag_presence[lag] = 1;
-          try
-            {
-              getDerivID(symbol_table.getID(eEndogenous, endoID), lag-1);
-            }
-          catch (UnknownDerivIDException &e)
-            {
-              lag_presence[lag] = 0;
-            }
-        }
-      if (lag_presence[0] == 1)
-        if (lag_presence[2] == 1)
-          zeta_mixed.push_back(endoID);
-        else
-          zeta_back.push_back(endoID);
-      else if (lag_presence[2] == 1)
-        zeta_fwrd.push_back(endoID);
-      else
-        zeta_static.push_back(endoID);
-
-    }
-  output << "size_t nstatic = " << zeta_static.size() << ";" << endl
-         << "size_t nfwrd   = " << zeta_fwrd.size() << ";" << endl
-         << "size_t nback   = " << zeta_back.size() << ";" << endl
-         << "size_t nmixed  = " << zeta_mixed.size() << ";" << endl;
-  output << "size_t zeta_static[" << zeta_static.size() << "] = {";
-  for (auto i = zeta_static.begin(); i != zeta_static.end(); ++i)
-    {
-      if (i != zeta_static.begin())
-        output << ",";
-      output << *i;
-    }
-  output << "};" << endl;
-
-  output << "size_t zeta_back[" << zeta_back.size() << "] = {";
-  for (auto i = zeta_back.begin(); i != zeta_back.end(); ++i)
-    {
-      if (i != zeta_back.begin())
-        output << ",";
-      output << *i;
-    }
-  output << "};" << endl;
-
-  output << "size_t zeta_fwrd[" << zeta_fwrd.size() << "] = {";
-  for (auto i = zeta_fwrd.begin(); i != zeta_fwrd.end(); ++i)
-    {
-      if (i != zeta_fwrd.begin())
-        output << ",";
-      output << *i;
-    }
-  output << "};" << endl;
-
-  output << "size_t zeta_mixed[" << zeta_mixed.size() << "] = {";
-  for (auto i = zeta_mixed.begin(); i != zeta_mixed.end(); ++i)
-    {
-      if (i != zeta_mixed.begin())
-        output << ",";
-      output << *i;
-    }
-  output << "};" << endl;
-
-  // Write number of non-zero derivatives
-  // Use -1 if the derivatives have not been computed
-  output << "int *NNZDerivatives[3] = {";
-  switch (order)
-    {
-    case 0:
-      output << NNZDerivatives[0] << ",-1,-1};" << endl;
-      break;
-    case 1:
-      output << NNZDerivatives[0] << "," << NNZDerivatives[1] << ",-1};" << endl;
-      break;
-    case 2:
-      output << NNZDerivatives[0] << "," << NNZDerivatives[1] << "," << NNZDerivatives[2] << "};" << endl;
-      break;
-    default:
-      cerr << "Order larger than 3 not implemented" << endl;
-      exit(EXIT_FAILURE);
-    }
-}
-
-void
-DynamicModel::writeResidualsC(const string &basename, bool cuda) const
-{
-  string filename = basename + "_residuals.c";
-  ofstream mDynamicModelFile, mDynamicMexFile;
-
-  mDynamicModelFile.open(filename, 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 residuals of the model for Dynare" << endl
-                    << " *" << endl
-                    << " * Warning : this file is generated automatically by Dynare" << endl
-                    << " *           from model " << basename << "(.mod)" << endl
-                    << " */" << endl
-#if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
-                    << "#ifdef _MSC_VER" << endl
-                    << "#define _USE_MATH_DEFINES" << endl
-                    << "#endif" << endl
-#endif
-                    << "#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
-  // even for residuals if doing Ramsey
-  writePowerDerivCHeader(mDynamicModelFile);
-  writeNormcdfCHeader(mDynamicModelFile);
-
-  mDynamicModelFile << "void Residuals(const double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual)" << endl
-                    << "{" << endl;
-
-  ostringstream model_output;    // Used for storing model equations
-  writeModelEquations(model_output, oCDynamic2Model);
-
-  mDynamicModelFile << "  double lhs, rhs;" << endl
-                    << endl
-                    << "  /* Residual equations */" << endl
-                    << model_output.str()
-                    << "}" << endl;
-
-  writePowerDeriv(mDynamicModelFile);
-  writeNormcdf(mDynamicModelFile);
-  mDynamicModelFile.close();
-
-}
-
-void
-DynamicModel::writeFirstDerivativesC(const string &basename, bool cuda) const
-{
-  string filename = basename + "_first_derivatives.c";
-  ofstream mDynamicModelFile, mDynamicMexFile;
-
-  mDynamicModelFile.open(filename, 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
-#if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
-                    << "#ifdef _MSC_VER" << endl
-                    << "#define _USE_MATH_DEFINES" << endl
-                    << "#endif" << endl
-#endif
-                    << "#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);
-  writeNormcdfCHeader(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;
-
-  // Writing Jacobian
-  for (const auto & first_derivative : first_derivatives)
-    {
-      int eq = first_derivative.first.first;
-      int var = first_derivative.first.second;
-      expr_t d1 = first_derivative.second;
-
-      jacobianHelper(mDynamicModelFile, eq, getDynJacobianCol(var), oCDynamicModel);
-      mDynamicModelFile << "=";
-      // oCStaticModel makes reference to the static variables
-      // oCDynamicModel makes reference to the dynamic variables
-      d1->writeOutput(mDynamicModelFile, oCDynamicModel, temporary_terms, {}, {});
-      mDynamicModelFile << ";" << endl;
-    }
-
-  mDynamicModelFile << "}" << endl;
-
-  mDynamicModelFile.close();
-
-}
-
-// using compressed sparse row format (CSR)
-void
-DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) const
-{
-  string filename = basename + "_first_derivatives.c";
-  ofstream mDynamicModelFile, mDynamicMexFile;
-
-  mDynamicModelFile.open(filename, 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
-#if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
-                    << "#ifdef _MSC_VER" << endl
-                    << "#define _USE_MATH_DEFINES" << endl
-                    << "#endif" << endl
-#endif
-                    << "#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);
-  writeNormcdfCHeader(mDynamicModelFile);
-
-  mDynamicModelFile << "void FirstDerivatives(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;
-
-  int cols_nbr = 3*symbol_table.endo_nbr() + symbol_table.exo_nbr() + symbol_table.exo_det_nbr();
-
-  // Indexing derivatives in column order
-  vector<derivative> D;
-  for (const auto & first_derivative : first_derivatives)
-    {
-      int eq = first_derivative.first.first;
-      int dynvar = first_derivative.first.second;
-      int lag = getLagByDerivID(dynvar);
-      int symb_id = getSymbIDByDerivID(dynvar);
-      SymbolType type = getTypeByDerivID(dynvar);
-      int tsid = symbol_table.getTypeSpecificID(symb_id);
-      int col_id;
-      switch (type)
-        {
-        case eEndogenous:
-          col_id = tsid+(lag+1)*symbol_table.endo_nbr();
-          break;
-        case eExogenous:
-          col_id = tsid+3*symbol_table.endo_nbr();
-          break;
-        case eExogenousDet:
-          col_id = tsid+3*symbol_table.endo_nbr()+symbol_table.exo_nbr();
-          break;
-        default:
-          std::cerr << "This case shouldn't happen" << std::endl;
-          exit(EXIT_FAILURE);
-        }
-      derivative deriv(col_id + eq *cols_nbr, col_id, eq, first_derivative.second);
-      D.push_back(deriv);
-    }
-  sort(D.begin(), D.end(), derivative_less_than());
-
-  // writing sparse Jacobian
-  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, oCDynamic2Model, temporary_terms, {}, {});
-      mDynamicModelFile << ";" << endl;
-      k++;
-    }
-
-  // row_ptr must point to the relative address of the first element of the row
-  int cumsum = 0;
-  mDynamicModelFile << "int row_ptr_data[" <<  row_ptr.size() + 1 << "] = { 0";
-  for (int & it : row_ptr)
-    {
-      cumsum += it;
-      mDynamicModelFile << ", " << cumsum;
-    }
-  mDynamicModelFile << "};" << endl
-                    << "int i;" << endl
-                    << "for (i=0; i < " << row_ptr.size() + 1 << "; i++) row_ptr[i] = row_ptr_data[i];" << endl;
-  mDynamicModelFile << "}" << endl;
-
-  mDynamicModelFile.close();
-
-}
-void
-DynamicModel::writeSecondDerivativesC_csr(const string &basename, bool cuda) const
-{
-
-  string filename = basename + "_second_derivatives.c";
-  ofstream mDynamicModelFile, mDynamicMexFile;
-
-  mDynamicModelFile.open(filename, 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
-#if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
-                    << "#ifdef _MSC_VER" << endl
-                    << "#define _USE_MATH_DEFINES" << endl
-                    << "#endif" << endl
-#endif
-                    << "#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);
-  writeNormcdfCHeader(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;
-
-  // Indexing derivatives in column order
-  vector<derivative> D;
-  int hessianColsNbr = dynJacobianColsNbr*dynJacobianColsNbr;
-  for (const auto & second_derivative : second_derivatives)
-    {
-      int eq, var1, var2;
-      tie(eq, var1, var2) = second_derivative.first;
-
-      int id1 = getDynJacobianCol(var1);
-      int id2 = getDynJacobianCol(var2);
-
-      int col_nb = id1 * dynJacobianColsNbr + id2;
-
-      derivative deriv(col_nb + eq *hessianColsNbr, col_nb, eq, second_derivative.second);
-      D.push_back(deriv);
-      if (id1 != id2)
-        {
-          col_nb = id2 * dynJacobianColsNbr + id1;
-          derivative deriv(col_nb + eq *hessianColsNbr, col_nb, eq, second_derivative.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, {}, {});
-      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 (int & it : row_ptr)
-    {
-      cumsum += it;
-      mDynamicModelFile << ", " << cumsum;
-    }
-  mDynamicModelFile << "];" << endl;
-
-  mDynamicModelFile << "}" << endl;
-
-  writePowerDeriv(mDynamicModelFile);
-  writeNormcdf(mDynamicModelFile);
-  mDynamicModelFile.close();
-}
-
-void
-DynamicModel::writeThirdDerivativesC_csr(const string &basename, bool cuda) const
-{
-  string filename = basename + "_third_derivatives.c";
-  ofstream mDynamicModelFile, mDynamicMexFile;
-
-  mDynamicModelFile.open(filename, 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
-#if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
-                    << "#ifdef _MSC_VER" << endl
-                    << "#define _USE_MATH_DEFINES" << endl
-                    << "#endif" << endl
-#endif
-                    << "#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);
-  writeNormcdfCHeader(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;
-
-  vector<derivative> D;
-  int hessianColsNbr = dynJacobianColsNbr*dynJacobianColsNbr;
-  int thirdDerivativesColsNbr = hessianColsNbr*dynJacobianColsNbr;
-  for (const auto & third_derivative : third_derivatives)
-    {
-      int eq, var1, var2, var3;
-      tie(eq, var1, var2, var3) = third_derivative.first;
-
-      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, third_derivative.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, third_derivative.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, third_derivative.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, third_derivative.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, third_derivative.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, third_derivative.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, {}, {});
-      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 (int & it : row_ptr)
-    {
-      cumsum += it;
-      mDynamicModelFile << ", " << cumsum;
-    }
-  mDynamicModelFile << "];" << endl;
-
-  mDynamicModelFile << "}" << endl;
-
-  writePowerDeriv(mDynamicModelFile);
-  writeNormcdf(mDynamicModelFile);
-  mDynamicModelFile.close();
-
-}
-
-void
-DynamicModel::writeCCOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present) const
-{
-  int lag_presence[3];
-  // Loop on endogenous variables
-  for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++)
-    {
-      // Loop on periods
-      for (int lag = 0; lag <= 2; lag++)
-        {
-          lag_presence[lag] = 1;
-          try
-            {
-              getDerivID(symbol_table.getID(eEndogenous, endoID), lag-1);
-            }
-          catch (UnknownDerivIDException &e)
-            {
-              lag_presence[lag] = 0;
-            }
-        }
-      if (lag_presence[0] == 1)
-        if (lag_presence[2] == 1)
-          output << "zeta_mixed.push_back(" << endoID << ");" << endl;
-        else
-          output << "zeta_back.push_back(" << endoID << ");" << endl;
-      else if (lag_presence[2] == 1)
-        output << "zeta_fwrd.push_back(" << endoID << ");" << endl;
-      else
-        output << "zeta_static.push_back(" << endoID << ");" << endl;
-
-    }
-  output << "nstatic = zeta_static.size();" << endl
-         << "nfwrd   = zeta_fwrd.size();" << endl
-         << "nback   = zeta_back.size();" << endl
-         << "nmixed  = zeta_mixed.size();" << endl;
-
-  // Write number of non-zero derivatives
-  // Use -1 if the derivatives have not been computed
-  output << endl
-         << "NNZDerivatives.push_back(" << NNZDerivatives[0] << ");" << endl;
-  if (order > 1)
-    {
-      output << "NNZDerivatives.push_back(" << NNZDerivatives[1] << ");" << endl;
-      if (order > 2)
-        output << "NNZDerivatives.push_back(" << NNZDerivatives[2] << ");" << endl;
-      else
-        output << "NNZDerivatives.push_back(-1);" << endl;
-    }
-  else
-    output << "NNZDerivatives.push_back(-1);" << endl
-           << "NNZDerivatives.push_back(-1);" << endl;
-}
-
 void
 DynamicModel::writeJsonOutput(ostream &output) const
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index cf38385ee7f8b7ebe556ca22e7e7df10a5c33010..80a93d2d9bcfd12a7979a9aad2458987f9556210 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -627,21 +627,6 @@ public:
   //! Returns true if a parameter was used in the model block with a lead or lag
   bool ParamUsedWithLeadLag() 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;
-  //! Writes C file containing residuals
-  void writeResidualsC(const string &basename, bool cuda) 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 first order derivatives of model evaluated at steady state (conpressed sparse column)
-  void writeFirstDerivativesC_csr(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;
-
   bool isChecksumMatching(const string &basename) const;
 };
 
diff --git a/src/DynareMain.cc b/src/DynareMain.cc
index 00974428e90c73f9d388890811d1d96324ec4628..974b331aaaea64d2c89b2f00fdc94fd1913fe3c5 100644
--- a/src/DynareMain.cc
+++ b/src/DynareMain.cc
@@ -57,7 +57,7 @@ usage()
 {
   cerr << "Dynare usage: dynare mod_file [debug] [noclearall] [onlyclearglobals] [savemacro[=macro_file]] [onlymacro] [nolinemacro] [noemptylinemacro] [notmpterms] [nolog] [warn_uninit]"
        << " [console] [nograph] [nointeractive] [parallel[=cluster_name]] [conffile=parallel_config_path_and_filename] [parallel_slave_open_mode] [parallel_test]"
-       << " [-D<variable>[=<value>]] [-I/path] [nostrict] [stochastic] [fast] [minimal_workspace] [compute_xrefs] [output=dynamic|first|second|third] [language=C|C++|julia]"
+       << " [-D<variable>[=<value>]] [-I/path] [nostrict] [stochastic] [fast] [minimal_workspace] [compute_xrefs] [output=dynamic|first|second|third] [language=julia]"
        << " [params_derivs_order=0|1|2] [transform_unary_ops]"
 #if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
        << " [cygwin] [msvc] [mingw]"
@@ -292,11 +292,7 @@ main(int argc, char **argv)
             {
               // we don't want temp terms in external functions (except Julia)
               no_tmp_terms = true;
-              if (strlen(argv[arg]) == 10 && !strncmp(argv[arg] + 9, "C", 1))
-                language = c;
-              else if (strlen(argv[arg]) ==  12 && !strncmp(argv[arg] + 9, "C++", 3))
-                language = cpp;
-              else if (strlen(argv[arg]) == 13 && !strncmp(argv[arg] + 9, "cuda", 4))
+              if (strlen(argv[arg]) == 13 && !strncmp(argv[arg] + 9, "cuda", 4))
                 language = cuda;
               else if (strlen(argv[arg]) == 15 && !strncmp(argv[arg] + 9, "python", 6))
                 language = python;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 4d542c99566d3f2ae4da345515d48d753544ed63..523668dcc5186404cacff7a6dfab9b1a60366e3a 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -864,10 +864,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           i = datatree.getDynJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type);
           output <<  "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
-        case oCDynamic2Model:
-          i = tsid + (lag+1)*datatree.symbol_table.endo_nbr() + ARRAY_SUBSCRIPT_OFFSET(output_type);
-          output <<  "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
-          break;
         case oCStaticModel:
         case oJuliaStaticModel:
         case oMatlabStaticModel:
@@ -899,9 +895,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case oSteadyStateFile:
           output << "ys_" << LEFT_ARRAY_SUBSCRIPT(output_type) << tsid + 1 << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
-        case oCSteadyStateFile:
-          output << "ys_[" << tsid << "]";
-          break;
         case oMatlabDseries:
           output << "ds." << datatree.symbol_table.getName(symb_id);
           if (lag != 0)
@@ -931,7 +924,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                    << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCDynamicModel:
-        case oCDynamic2Model:
           if (lag == 0)
             output <<  "x[it_+" << i << "*nb_row_x]";
           else if (lag > 0)
@@ -956,9 +948,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case oSteadyStateFile:
           output << "exo_" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
-        case oCSteadyStateFile:
-          output << "exo_[" << i - 1 << "]";
-          break;
         case oMatlabDseries:
           output << "ds." << datatree.symbol_table.getName(symb_id);
           if (lag != 0)
@@ -988,7 +977,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
                    << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCDynamicModel:
-        case oCDynamic2Model:
           if (lag == 0)
             output <<  "x[it_+" << i << "*nb_row_x]";
           else if (lag > 0)
@@ -1013,9 +1001,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case oSteadyStateFile:
           output << "exo_" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
-        case oCSteadyStateFile:
-          output << "exo_[" << i - 1 << "]";
-          break;
         case oMatlabDseries:
           output << "ds." << datatree.symbol_table.getName(symb_id);
           if (lag != 0)
@@ -6686,7 +6671,7 @@ ExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_typ
                                   const deriv_node_temp_terms_t &tef_terms) const
 {
   if (output_type == oMatlabOutsideModel || output_type == oSteadyStateFile
-      || output_type == oCSteadyStateFile || output_type == oJuliaSteadyStateFile
+      || output_type == oJuliaSteadyStateFile
       || IS_LATEX(output_type))
     {
       string name = IS_LATEX(output_type) ? datatree.symbol_table.getTeXName(symb_id)
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 44999bc22c1dc00dd5203898dacd3e0ea7deb536..e14f2ad2d4df1ae03622c34d0686fb903a26888b 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -74,7 +74,6 @@ enum ExprNodeOutputType
     oMatlabStaticModelSparse,                     //!< Matlab code, static block decomposed model
     oMatlabDynamicModelSparse,                    //!< Matlab code, dynamic block decomposed model
     oCDynamicModel,                               //!< C code, dynamic model
-    oCDynamic2Model,                              //!< C code, dynamic model, alternative numbering of endogenous variables
     oCStaticModel,                                //!< C code, static model
     oJuliaStaticModel,                            //!< Julia code, static model
     oJuliaDynamicModel,                           //!< Julia code, dynamic model
@@ -87,7 +86,6 @@ enum ExprNodeOutputType
     oCDynamicSteadyStateOperator,                 //!< C code, dynamic model, inside a steady state operator
     oJuliaDynamicSteadyStateOperator,             //!< Julia code, dynamic model, inside a steady state operator
     oSteadyStateFile,                             //!< Matlab code, in the generated steady state file
-    oCSteadyStateFile,                            //!< C code, in the generated steady state file
     oJuliaSteadyStateFile,                        //!< Julia code, in the generated steady state file
     oMatlabDseries                                //!< Matlab code for dseries
   };
@@ -108,10 +106,8 @@ enum ExprNodeOutputType
                                || (output_type) == oJuliaSteadyStateFile)
 
 #define IS_C(output_type) ((output_type) == oCDynamicModel              \
-                           || (output_type) == oCDynamic2Model          \
                            || (output_type) == oCStaticModel            \
-                           || (output_type) == oCDynamicSteadyStateOperator \
-                           || (output_type) == oCSteadyStateFile)
+                           || (output_type) == oCDynamicSteadyStateOperator)
 
 #define IS_LATEX(output_type) ((output_type) == oLatexStaticModel       \
                                || (output_type) == oLatexDynamicModel   \
diff --git a/src/ExtendedPreprocessorTypes.hh b/src/ExtendedPreprocessorTypes.hh
index e40c73eeb5f7794f1e7daf94eee0a05e4b5b40b9..87fea10cf94035463bdc9b608c70515ea7f1e25d 100644
--- a/src/ExtendedPreprocessorTypes.hh
+++ b/src/ExtendedPreprocessorTypes.hh
@@ -32,8 +32,6 @@ enum FileOutputType
 enum LanguageOutputType
   {
     matlab,                           // outputs files for Matlab/Octave processing
-    c,                                // outputs files for C
-    cpp,                              // outputs files for C++
     cuda,                             // outputs files for CUDA (not yet implemented)
     julia,                            // outputs files for Julia
     python,                           // outputs files for Python (not yet implemented) (not yet implemented)
diff --git a/src/ModFile.cc b/src/ModFile.cc
index c0d932d89f5aac1ee484d6ef7954a506c234c33b..62f892f78e8cb89b1b3bac8d5bddb1369c990483 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -1061,12 +1061,6 @@ ModFile::writeExternalFiles(const string &basename, FileOutputType output, Langu
 {
   switch (language)
     {
-    case c:
-      writeExternalFilesC(basename, output);
-      break;
-    case cpp:
-      writeExternalFilesCC(basename, output);
-      break;
     case julia:
       writeExternalFilesJulia(basename, output, nopreprocessoroutput);
       break;
@@ -1076,213 +1070,6 @@ ModFile::writeExternalFiles(const string &basename, FileOutputType output, Langu
     }
 }
 
-void
-ModFile::writeExternalFilesC(const string &basename, FileOutputType output) const
-{
-  writeModelC(basename);
-  steady_state_model.writeSteadyStateFileC(basename, mod_file_struct.ramsey_model_present);
-
-  dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option, false);
-
-  if (!no_static)
-    static_model.writeStaticFile(basename, false, false, true, false);
-
-  //  static_model.writeStaticCFile(basename, block, byte_code, use_dll);
-  //  static_model.writeParamsDerivativesFileC(basename, cuda);
-  //  static_model.writeAuxVarInitvalC(mOutputFile, oMatlabOutsideModel, cuda);
-
-  dynamic_model.writeResidualsC(basename, cuda);
-  // dynamic_model.writeParamsDerivativesFileC(basename, cuda);
-  dynamic_model.writeFirstDerivativesC(basename, cuda);
-
-  if (output == second)
-    dynamic_model.writeSecondDerivativesC_csr(basename, cuda);
-  else if (output == third)
-    {
-      dynamic_model.writeSecondDerivativesC_csr(basename, cuda);
-      dynamic_model.writeThirdDerivativesC_csr(basename, cuda);
-    }
-}
-
-void
-ModFile::writeModelC(const string &basename) const
-{
-  string filename = basename + ".c";
-
-  ofstream mDriverCFile;
-  mDriverCFile.open(filename, ios::out | ios::binary);
-  if (!mDriverCFile.is_open())
-    {
-      cerr << "Error: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  mDriverCFile << "/*" << endl
-               << " * " << filename << " : Driver file for Dynare C code" << endl
-               << " *" << endl
-               << " * Warning : this file is generated automatically by Dynare" << endl
-               << " *           from model file (.mod)" << endl
-               << " */" << endl
-               << endl
-               << "#include \"dynare_driver.h\"" << endl
-               << endl
-               << "struct" << endl
-               << "{" << endl;
-
-  // Write basic info
-  symbol_table.writeCOutput(mDriverCFile);
-
-  mDriverCFile << endl << "params.resize(param_nbr);" << endl;
-
-  if (dynamic_model.equation_number() > 0)
-    {
-      dynamic_model.writeCOutput(mDriverCFile, basename, block, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present);
-      //      if (!no_static)
-      //        static_model.writeCOutput(mOutputFile, block);
-    }
-
-  // Print statements
-  for (auto statement : statements)
-    statement->writeCOutput(mDriverCFile, basename);
-
-  mDriverCFile << "} DynareInfo;" << endl;
-  mDriverCFile.close();
-
-  // Write informational m file
-  ofstream mOutputFile;
-
-  if (basename.size())
-    {
-      string fname(basename);
-      fname += ".m";
-      mOutputFile.open(fname, ios::out | ios::binary);
-      if (!mOutputFile.is_open())
-        {
-          cerr << "ERROR: Can't open file " << fname
-               << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-    }
-  else
-    {
-      cerr << "ERROR: Missing file name" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  mOutputFile << "%" << endl
-              << "% Status : informational m file" << endl
-              << "%" << endl
-              << "% Warning : this file is generated automatically by Dynare" << endl
-              << "%           from model file (.mod)" << endl << endl
-              << "disp('The following C file was successfully created:');" << endl
-              << "ls preprocessorOutput.c" << endl << endl;
-  mOutputFile.close();
-}
-
-void
-ModFile::writeExternalFilesCC(const string &basename, FileOutputType output) const
-{
-  writeModelCC(basename);
-  steady_state_model.writeSteadyStateFileC(basename, mod_file_struct.ramsey_model_present);
-
-  dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option, false);
-
-  if (!no_static)
-    static_model.writeStaticFile(basename, false, false, true, false);
-
-  //  static_model.writeStaticCFile(basename, block, byte_code, use_dll);
-  //  static_model.writeParamsDerivativesFileC(basename, cuda);
-  //  static_model.writeAuxVarInitvalC(mOutputFile, oMatlabOutsideModel, cuda);
-
-  // dynamic_model.writeResidualsC(basename, cuda);
-  // dynamic_model.writeParamsDerivativesFileC(basename, cuda);
-  dynamic_model.writeResidualsC(basename, cuda);
-  dynamic_model.writeFirstDerivativesC_csr(basename, cuda);
-
-  if (output == second)
-    dynamic_model.writeSecondDerivativesC_csr(basename, cuda);
-  else if (output == third)
-    {
-      dynamic_model.writeSecondDerivativesC_csr(basename, cuda);
-      dynamic_model.writeThirdDerivativesC_csr(basename, cuda);
-    }
-}
-
-void
-ModFile::writeModelCC(const string &basename) const
-{
-  string filename = basename + ".cc";
-
-  ofstream mDriverCFile;
-  mDriverCFile.open(filename, ios::out | ios::binary);
-  if (!mDriverCFile.is_open())
-    {
-      cerr << "Error: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  mDriverCFile << "/*" << endl
-               << " * " << filename << " : Driver file for Dynare C++ code" << endl
-               << " *" << endl
-               << " * Warning : this file is generated automatically by Dynare" << endl
-               << " *           from model file (.mod)" << endl
-               << " */" << endl
-               << endl
-               << "#include \"dynare_cpp_driver.hh\"" << endl
-               << endl
-               << "DynareInfo::DynareInfo(void)" << endl
-               << "{" << endl;
-
-  // Write basic info
-  symbol_table.writeCCOutput(mDriverCFile);
-
-  mDriverCFile << endl << "params.resize(param_nbr);" << endl;
-
-  if (dynamic_model.equation_number() > 0)
-    {
-      dynamic_model.writeCCOutput(mDriverCFile, basename, block, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present);
-      //      if (!no_static)
-      //        static_model.writeCOutput(mOutputFile, block);
-    }
-
-  // Print statements
-  for (auto statement : statements)
-    statement->writeCOutput(mDriverCFile, basename);
-
-  mDriverCFile << "};" << endl;
-  mDriverCFile.close();
-
-  // Write informational m file
-  ofstream mOutputFile;
-
-  if (basename.size())
-    {
-      string fname(basename);
-      fname += ".m";
-      mOutputFile.open(fname, ios::out | ios::binary);
-      if (!mOutputFile.is_open())
-        {
-          cerr << "ERROR: Can't open file " << fname
-               << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-    }
-  else
-    {
-      cerr << "ERROR: Missing file name" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  mOutputFile << "%" << endl
-              << "% Status : informational m file" << endl
-              << "%" << endl
-              << "% Warning : this file is generated automatically by Dynare" << endl
-              << "%           from model file (.mod)" << endl << endl
-              << "disp('The following C++ file was successfully created:');" << endl
-              << "ls preprocessorOutput.cc" << endl << endl;
-  mOutputFile.close();
-}
-
 void
 ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output, const bool nopreprocessoroutput) const
 {
diff --git a/src/ModFile.hh b/src/ModFile.hh
index d22203cbce989bdcc8bc4ac8327186ffc7ba263a..1ecac964fceeb457797667f3821f4eec22cabaf8 100644
--- a/src/ModFile.hh
+++ b/src/ModFile.hh
@@ -157,15 +157,7 @@ public:
                         , const bool nopreprocessoroutput
                         ) const;
   void writeExternalFiles(const string &basename, FileOutputType output, LanguageOutputType language, const bool nopreprocessoroutput) const;
-  void writeExternalFilesC(const string &basename, FileOutputType output) const;
-  void writeExternalFilesCC(const string &basename, FileOutputType output) const;
   void writeExternalFilesJulia(const string &basename, FileOutputType output, const bool nopreprocessoroutput) const;
-  //! Writes C output files only => No further Matlab processing
-  void writeCOutputFiles(const string &basename) const;
-  void writeModelC(const string &basename) const;
-  //! Writes Cpp output files only => No further Matlab processing
-  void writeCCOutputFiles(const string &basename) const;
-  void writeModelCC(const string &basename) const;
 
   void computeChecksum();
   //! Write JSON representation of ModFile object
diff --git a/src/SteadyStateModel.cc b/src/SteadyStateModel.cc
index 326a64c1016438fee748c120b1cf0515b84fc2d0..215b140eaafdd8be1f373751c7f16454ebb25407 100644
--- a/src/SteadyStateModel.cc
+++ b/src/SteadyStateModel.cc
@@ -220,54 +220,6 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
     output << "end" << endl;
 }
 
-void
-SteadyStateModel::writeSteadyStateFileC(const string &basename, bool ramsey_model) const
-{
-  string filename = basename + "_steadystate.c";
-
-  ofstream output;
-  output.open(filename, 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 (const auto & i : def_table)
-    {
-      const vector<int> &symb_ids = i.first;
-      output << "    ";
-      if (symb_ids.size() > 1)
-        std::cout << "Error: in C, multiple returns are not permitted in steady_state_model" << std::endl;
-      auto it = variable_node_map.find({ 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 << "=";
-      i.second->writeOutput(output, oCSteadyStateFile);
-      output << ";" << endl;
-    }
-  output << "    // Auxiliary equations" << endl;
-  static_model.writeAuxVarInitval(output, oCSteadyStateFile);
-  output << "}" << endl;
-}
-
 void
 SteadyStateModel::writeJsonSteadyStateFile(ostream &output, bool transformComputingPass) const
 {
diff --git a/src/SteadyStateModel.hh b/src/SteadyStateModel.hh
index 977765049e9aa2b4ba3dba0c8192df1870e57246..04166dd80464eca09bf837db5648c4872b37899a 100644
--- a/src/SteadyStateModel.hh
+++ b/src/SteadyStateModel.hh
@@ -50,7 +50,6 @@ 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, bool julia) const;
-  void writeSteadyStateFileC(const string &basename, bool ramsey_model) const;
   //! Writes LaTeX file with the equations of the dynamic model (for the steady state model)
   void writeLatexSteadyStateFile(const string &basename) const;
   //! Writes JSON output