diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index f21c6264826d21ece3813684739d1e8c596dd976..55f804cfdb04b479a12ea6995ab9344c714f1a1b 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -950,140 +950,6 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
   writeToFileIfModified(output, basename + "Dynamic.jl");
 }
 
-void
-DynamicModel::writeDynamicCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const
-{
-  string filename = basename + "/model/src/dynamic.c";
-
-  int ntt { static_cast<int>(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size()) };
-
-  ofstream output{filename, ios::out | ios::binary};
-  if (!output.is_open())
-    {
-      cerr << "Error: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-  output << "/*" << endl
-         << " * " << filename << " : Computes " << modelClassName() << " for Dynare" << endl
-         << " *" << endl
-         << " * Warning : this file is generated automatically by Dynare" << endl
-         << " *           from model file (.mod)" << endl
-         << " */" << endl
-         << endl
-         << "#include <math.h>" << endl
-         << "#include <stdlib.h>" << endl
-         << R"(#include "mex.h")" << endl
-         << endl;
-
-  // Write function definition if BinaryOpcode::powerDeriv is used
-  writePowerDeriv(output);
-
-  output << endl;
-
-  auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::CDynamicModel>();
-
-  for (size_t i{0}; i < d_output.size(); i++)
-    {
-      string funcname{i == 0 ? "resid" : "g" + to_string(i)};
-      output << "void dynamic_" << funcname << "_tt(const double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, int it_, double *restrict T)" << endl
-             << "{" << endl
-             << tt_output[i].str()
-             << "}" << endl
-             << endl
-             << "void dynamic_" << funcname << "(const double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, int it_, const double *restrict T, ";
-      if (i == 0)
-        output << "double *restrict residual";
-      else if (i == 1)
-        output << "double *restrict g1";
-      else
-        output << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v";
-      output << ")" << endl
-             << "{" << endl;
-      if (i == 0)
-        output << "  double lhs, rhs;" << endl;
-      output << d_output[i].str()
-             << "}" << endl
-             << endl;
-    }
-
-  output << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
-         << "{" << endl
-         << "  if (nlhs > " << min(computed_derivs_order + 1, 4) << ")" << endl
-         << R"(    mexErrMsgTxt("Derivatives of higher order than computed have been requested");)" << endl
-         << "  if (nrhs != 5)" << endl
-         << R"(    mexErrMsgTxt("Requires exactly 5 input arguments");)" << endl
-         << endl
-         << "  double *y = mxGetPr(prhs[0]);" << endl
-         << "  double *x = mxGetPr(prhs[1]);" << endl
-         << "  double *params = mxGetPr(prhs[2]);" << endl
-         << "  double *steady_state = mxGetPr(prhs[3]);" << endl
-         << "  int it_ = (int) mxGetScalar(prhs[4]) - 1;" << endl
-         << "  int nb_row_x = mxGetM(prhs[1]);" << endl
-         << endl
-         << "  double *T = (double *) malloc(sizeof(double)*" << ntt << ");" << endl
-         << endl
-         << "  if (nlhs >= 1)" << endl
-         << "    {" << endl
-         << "       plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
-         << "       double *residual = mxGetPr(plhs[0]);" << endl
-         << "       dynamic_resid_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl
-         << "       dynamic_resid(y, x, nb_row_x, params, steady_state, it_, T, residual);" << endl
-         << "    }" << endl
-         << endl
-         << "  if (nlhs >= 2)" << endl
-         << "    {" << endl
-         << "       plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << getJacobianColsNbr() << ", mxREAL);" << endl
-         << "       double *g1 = mxGetPr(plhs[1]);" << endl
-         << "       dynamic_g1_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl
-         << "       dynamic_g1(y, x, nb_row_x, params, steady_state, it_, T, g1);" << endl
-         << "    }" << endl
-         << endl
-         << "  if (nlhs >= 3)" << endl
-         << "    {" << endl
-         << "      mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      dynamic_g2_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl
-         << "      dynamic_g2(y, x, nb_row_x, params, steady_state, it_, T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
-         << "      mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
-         << "      mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr()*getJacobianColsNbr() << ");" << endl
-         << "      mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
-         << R"(      mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
-         << "      plhs[2] = plhs_sparse[0];" << endl
-         << "      mxDestroyArray(g2_i);" << endl
-         << "      mxDestroyArray(g2_j);" << endl
-         << "      mxDestroyArray(g2_v);" << endl
-         << "      mxDestroyArray(m);" << endl
-         << "      mxDestroyArray(n);" << endl
-         << "    }" << endl
-         << endl
-         << "  if (nlhs >= 4)" << endl
-         << "    {" << endl
-         << "      mxArray *g3_i = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
-         << "      mxArray *g3_j = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
-         << "      mxArray *g3_v = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
-         << "      dynamic_g3_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl
-         << "      dynamic_g3(y, x, nb_row_x, params, steady_state, it_, T, mxGetPr(g3_i), mxGetPr(g3_j), mxGetPr(g3_v));" << endl
-         << "      mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
-         << "      mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr()*getJacobianColsNbr()*getJacobianColsNbr() << ");" << endl
-         << "      mxArray *plhs_sparse[1], *prhs_sparse[5] = { g3_i, g3_j, g3_v, m, n };" << endl
-         << R"(      mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
-         << "      plhs[3] = plhs_sparse[0];" << endl
-         << "      mxDestroyArray(g3_i);" << endl
-         << "      mxDestroyArray(g3_j);" << endl
-         << "      mxDestroyArray(g3_v);" << endl
-         << "      mxDestroyArray(m);" << endl
-         << "      mxDestroyArray(n);" << endl
-         << "    }" << endl
-         << endl
-         << "  free(T);" << endl
-         << "}" << endl;
-
-  output.close();
-
-  compileMEX("+" + basename, "dynamic", mexext, { filename }, matlabroot, dynareroot);
-}
-
 string
 DynamicModel::reform(const string &name1) const
 {
@@ -3633,7 +3499,7 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool use_dll,
       writeDynamicBytecode(basename);
 
       if (use_dll)
-        writeDynamicCFile(basename, mexext, matlabroot, dynareroot);
+        writeModelCFile<true>(basename, mexext, matlabroot, dynareroot);
       else if (julia)
         writeDynamicJuliaFile(basename);
       else
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 1f0575889ffc08f0f559cba582b6dff146a184a0..7887c82fdbc7f7f4ea26080255c721f593838e75 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -121,8 +121,6 @@ private:
   void writeDynamicMFile(const string &basename) const;
   //! Writes dynamic model file (Julia version)
   void writeDynamicJuliaFile(const string &basename) const;
-  // Writes and compiles dynamic model file (C version)
-  void writeDynamicCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
   //! Writes the main dynamic function of block decomposed model (MATLAB version)
   void writeDynamicBlockMFile(const string &basename) const;
   /* Writes the main dynamic functions of block decomposed model (C version),
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 540209d48ea8686ca188d4b62a36aa6b504a82b2..08bf4d81c50b6d46bb56c0ff2afebc7bd08024ee 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -273,6 +273,10 @@ protected:
   template<ExprNodeOutputType output_type>
   pair<vector<ostringstream>, vector<ostringstream>> writeModelFileHelper() const;
 
+  // Writes and compiles dynamic/static file (C version)
+  template<bool dynamic>
+  void writeModelCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
+
   // Writes per-block residuals and temporary terms (incl. for derivatives)
   template<ExprNodeOutputType output_type>
   void writePerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
@@ -821,6 +825,216 @@ ModelTree::writeModelFileHelper() const
   return { move(d_output), move(tt_output) };
 }
 
+template<bool dynamic>
+void
+ModelTree::writeModelCFile(const string &basename, const string &mexext,
+                           const filesystem::path &matlabroot,
+                           const filesystem::path &dynareroot) const
+{
+  ofstream output;
+  auto open_file = [&output](const filesystem::path &p)
+  {
+    output.open(p, ios::out | ios::binary);
+    if (!output.is_open())
+      {
+        cerr << "ERROR: Can't open file " << p << " for writing" << endl;
+        exit(EXIT_FAILURE);
+      }
+  };
+
+  const filesystem::path model_src_dir { filesystem::path{basename} / "model" / "src" };
+
+  auto [d_output, tt_output] = writeModelFileHelper<dynamic ? ExprNodeOutputType::CDynamicModel : ExprNodeOutputType::CStaticModel>();
+  vector<filesystem::path> header_files, object_files;
+
+  // TODO: when C++20 support is complete, mark the following strings constexpr
+  const string prefix { dynamic ? "dynamic_" : "static_" };
+  const string ss_it_argin { dynamic ? ", const double *restrict steady_state, int it_" : "" };
+  const string ss_it_argout { dynamic ? ", steady_state, it_" : "" };
+  const string nb_row_x_argin { dynamic ? ", int nb_row_x" : "" };
+  const string nb_row_x_argout { dynamic ? ", nb_row_x" : "" };
+
+  for (size_t i {0}; i < d_output.size(); i++)
+    {
+      const string funcname { prefix + (i == 0 ? "resid" : "g" + to_string(i))};
+
+      const string prototype_tt { "void " + funcname + "_tt(const double *restrict y, const double *restrict x" + nb_row_x_argin + ", const double *restrict params" + ss_it_argin + ", double *restrict T)" };
+
+      const filesystem::path header_tt { model_src_dir / (funcname + "_tt.h") };
+      open_file(header_tt);
+      output << prototype_tt << ";" << endl;
+      output.close();
+      header_files.push_back(header_tt);
+
+      const filesystem::path source_tt { model_src_dir / (funcname + "_tt.c") };
+      open_file(source_tt);
+      output << "#include <math.h>" << endl
+             << R"(#include "mex.h")" << endl // Needed for calls to external functions
+             << endl;
+      writePowerDerivHeader(output);
+      output << endl
+             << prototype_tt << endl
+             << "{" << endl
+             << tt_output[i].str()
+             << "}" << endl
+             << endl;
+      output.close();
+      object_files.push_back(compileMEX(model_src_dir, funcname + "_tt" , mexext, { source_tt },
+                                        matlabroot, dynareroot, false));
+
+      const string prototype_main
+        {
+          [&funcname, &ss_it_argin, &nb_row_x_argin, i]
+          {
+            string p = "void " + funcname + "(const double *restrict y, const double *restrict x" + nb_row_x_argin + ", const double *restrict params" + ss_it_argin + ", const double *restrict T, ";
+            if (i == 0)
+              p += "double *restrict residual";
+            else if (i == 1)
+              p += "double *restrict g1";
+            else
+              p += "double *restrict g" + to_string(i) + "_i, double *restrict g" +
+                to_string(i) + "_j, double *restrict g" + to_string(i) + "_v";
+            p += ")";
+            return p;
+          }()
+        };
+
+      const filesystem::path header_main { model_src_dir / (funcname + ".h") };
+      open_file(header_main);
+      output << prototype_main << ";" << endl;
+      output.close();
+      header_files.push_back(header_main);
+
+      const filesystem::path source_main { model_src_dir / (funcname + ".c") };
+      open_file(source_main);
+      output << "#include <math.h>" << endl
+             << R"(#include "mex.h")" << endl // Needed for calls to external functions
+             << endl;
+      writePowerDerivHeader(output);
+      output << endl
+             << prototype_main << endl
+             << "{" << endl;
+      if (i == 0)
+        output << "  double lhs, rhs;" << endl;
+      output << d_output[i].str()
+             << "}" << endl
+             << endl;
+      output.close();
+      object_files.push_back(compileMEX(model_src_dir, funcname, mexext, { source_main },
+                                        matlabroot, dynareroot, false));
+    }
+
+  const filesystem::path filename { model_src_dir / (dynamic ? "dynamic.c" : "static.c") };
+
+  const int ntt { static_cast<int>(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size()) };
+
+  open_file(filename);
+  output << "/*" << endl
+         << " * " << filename << " : Computes " << modelClassName() << " for Dynare" << endl
+         << " *" << endl
+         << " * Warning : this file is generated automatically by Dynare" << endl
+         << " *           from model file (.mod)" << endl
+         << " */" << endl
+         << endl
+         << "#include <math.h>" << endl // Needed for getPowerDeriv()
+         << "#include <stdlib.h>" << endl // Needed for malloc() and free()
+         << R"(#include "mex.h")" << endl;
+  for (const auto &it : header_files)
+    output << "#include " << it.filename() << endl;
+  output << endl;
+
+  // Write function definition if BinaryOpcode::powerDeriv is used
+  writePowerDeriv(output);
+
+  output << endl
+         << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
+         << "{" << endl;
+  if constexpr(dynamic)
+    output << "  if (nlhs > " << min(computed_derivs_order + 1, 4) << ")" << endl
+           << R"(    mexErrMsgTxt("Derivatives of higher order than computed have been requested");)" << endl
+           << "  if (nrhs != 5)" << endl
+           << R"(    mexErrMsgTxt("Requires exactly 5 input arguments");)" << endl;
+  else
+    output << "  if (nrhs > 3)" << endl
+           << R"(    mexErrMsgTxt("Accepts at most 3 output arguments");)" << endl
+           << "  if (nrhs != 3)" << endl
+           << R"(    mexErrMsgTxt("Requires exactly 3 input arguments");)" << endl;
+  output << endl
+         << "  double *y = mxGetPr(prhs[0]);" << endl
+         << "  double *x = mxGetPr(prhs[1]);" << endl
+         << "  double *params = mxGetPr(prhs[2]);" << endl;
+  if constexpr(dynamic)
+    output << "  double *steady_state = mxGetPr(prhs[3]);" << endl
+           << "  int it_ = (int) mxGetScalar(prhs[4]) - 1;" << endl
+           << "  int nb_row_x = mxGetM(prhs[1]);" << endl;
+  output << endl
+         << "  double *T = (double *) malloc(sizeof(double)*" << ntt << ");" << endl
+         << endl
+         << "  if (nlhs >= 1)" << endl
+         << "    {" << endl
+         << "       plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
+         << "       double *residual = mxGetPr(plhs[0]);" << endl
+         << "       " << prefix << "resid_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T);" << endl
+         << "       " << prefix << "resid(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T, residual);" << endl
+         << "    }" << endl
+         << endl
+         << "  if (nlhs >= 2)" << endl
+         << "    {" << endl
+         << "       plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << getJacobianColsNbr() << ", mxREAL);" << endl
+         << "       double *g1 = mxGetPr(plhs[1]);" << endl
+         << "       " << prefix << "g1_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T);" << endl
+         << "       " << prefix << "g1(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T, g1);" << endl
+         << "    }" << endl
+         << endl
+         << "  if (nlhs >= 3)" << endl
+         << "    {" << endl
+         << "      mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+         << "      mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+         << "      mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
+         << "      " << prefix << "g2_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T);" << endl
+         << "      " << prefix << "g2(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
+         << "      mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
+         << "      mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr()*getJacobianColsNbr() << ");" << endl
+         << "      mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
+         << R"(      mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
+         << "      plhs[2] = plhs_sparse[0];" << endl
+         << "      mxDestroyArray(g2_i);" << endl
+         << "      mxDestroyArray(g2_j);" << endl
+         << "      mxDestroyArray(g2_v);" << endl
+         << "      mxDestroyArray(m);" << endl
+         << "      mxDestroyArray(n);" << endl
+         << "    }" << endl
+         << endl;
+  if constexpr(dynamic)
+    output << "  if (nlhs >= 4)" << endl
+           << "    {" << endl
+           << "      mxArray *g3_i = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
+           << "      mxArray *g3_j = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
+           << "      mxArray *g3_v = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1 << ", mxREAL);" << endl
+           << "      " << prefix << "g3_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T);" << endl
+           << "      " << prefix << "g3(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T, mxGetPr(g3_i), mxGetPr(g3_j), mxGetPr(g3_v));" << endl
+           << "      mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
+           << "      mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr()*getJacobianColsNbr()*getJacobianColsNbr() << ");" << endl
+           << "      mxArray *plhs_sparse[1], *prhs_sparse[5] = { g3_i, g3_j, g3_v, m, n };" << endl
+           << R"(      mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
+           << "      plhs[3] = plhs_sparse[0];" << endl
+           << "      mxDestroyArray(g3_i);" << endl
+           << "      mxDestroyArray(g3_j);" << endl
+           << "      mxDestroyArray(g3_v);" << endl
+           << "      mxDestroyArray(m);" << endl
+           << "      mxDestroyArray(n);" << endl
+           << "    }" << endl
+           << endl;
+
+  output << "  free(T);" << endl
+         << "}" << endl;
+  output.close();
+
+  object_files.push_back(filename);
+  compileMEX("+" + basename, dynamic ? "dynamic" : "static", mexext, object_files, matlabroot,
+             dynareroot);
+}
+
 template<ExprNodeOutputType output_type>
 void
 ModelTree::writePerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index ce63069c995dc248a83368379d8944a187e9c179..dde18cb9ccedd53f6caff370bceff911c4d832c6 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -610,119 +610,6 @@ StaticModel::writeStaticMCompatFile(const string &basename) const
   output.close();
 }
 
-void
-StaticModel::writeStaticCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const
-{
-  // Writing comments and function definition command
-  string filename{basename + "/model/src/static.c"};
-
-  int ntt { static_cast<int>(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size()) };
-
-  ofstream output{filename, ios::out | ios::binary};
-  if (!output.is_open())
-    {
-      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  output << "/*" << endl
-         << " * " << filename << " : Computes " << modelClassName() << " for Dynare" << endl
-         << " *" << endl
-         << " * Warning : this file is generated automatically by Dynare" << endl
-         << " *           from model file (.mod)" << endl << endl
-         << " */" << endl
-         << endl
-         << "#include <math.h>" << endl
-         << "#include <stdlib.h>" << endl
-         << R"(#include "mex.h")" << endl
-         << endl;
-
-  // Write function definition if BinaryOpcode::powerDeriv is used
-  writePowerDeriv(output);
-
-  output << endl;
-
-  auto [d_output, tt_output] = writeModelFileHelper<ExprNodeOutputType::CStaticModel>();
-
-  for (size_t i = 0; i < d_output.size(); i++)
-    {
-      string funcname{i == 0 ? "resid" : "g" + to_string(i)};
-      output << "void static_" << funcname << "_tt(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl
-             << "{" << endl
-             << tt_output[i].str()
-             << "}" << endl
-             << endl
-             << "void static_" << funcname << "(const double *restrict y, const double *restrict x, const double *restrict params, const double *restrict T, ";
-      if (i == 0)
-        output << "double *restrict residual";
-      else if (i == 1)
-        output << "double *restrict g1";
-      else
-        output << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v";
-      output << ")" << endl
-                   << "{" << endl;
-      if (i == 0)
-        output << "  double lhs, rhs;" << endl;
-      output << d_output[i].str()
-             << "}" << endl
-             << endl;
-    }
-
-  output << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
-         << "{" << endl
-         << "  if (nrhs > 3)" << endl
-         << R"(    mexErrMsgTxt("Accepts at most 3 output arguments");)" << endl
-         << "  if (nrhs != 3)" << endl
-         << R"(    mexErrMsgTxt("Requires exactly 3 input arguments");)" << endl
-         << "  double *y = mxGetPr(prhs[0]);" << endl
-         << "  double *x = mxGetPr(prhs[1]);" << endl
-         << "  double *params = mxGetPr(prhs[2]);" << endl
-         << endl
-         << "  double *T = (double *) malloc(sizeof(double)*" << ntt << ");" << endl
-         << endl
-         << "  if (nlhs >= 1)" << endl
-         << "    {" << endl
-         << "      plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
-         << "      double *residual = mxGetPr(plhs[0]);" << endl
-         << "      static_resid_tt(y, x, params, T);" << endl
-         << "      static_resid(y, x, params, T, residual);" << endl
-         << "    }" << endl
-         << endl
-         << "  if (nlhs >= 2)" << endl
-         << "    {" << endl
-         << "      plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << symbol_table.endo_nbr() << ", mxREAL);" << endl
-         << "      double *g1 = mxGetPr(plhs[1]);" << endl
-         << "      static_g1_tt(y, x, params, T);" << endl
-         << "      static_g1(y, x, params, T, g1);" << endl
-         << "    }" << endl
-         << endl
-         << "  if (nlhs >= 3)" << endl
-         << "    {" << endl
-         << "      mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1 << ", mxREAL);" << endl
-         << "      static_g2_tt(y, x, params, T);" << endl
-         << "      static_g2(y, x, params, T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
-         << "      mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
-         << "      mxArray *n = mxCreateDoubleScalar(" << symbol_table.endo_nbr()*symbol_table.endo_nbr() << ");" << endl
-         << "      mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
-         << R"(      mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
-         << "      plhs[2] = plhs_sparse[0];" << endl
-         << "      mxDestroyArray(g2_i);" << endl
-         << "      mxDestroyArray(g2_j);" << endl
-         << "      mxDestroyArray(g2_v);" << endl
-         << "      mxDestroyArray(m);" << endl
-         << "      mxDestroyArray(n);" << endl
-         << "    }" << endl
-         << endl
-         << "  free(T);" << endl
-         << "}" << endl;
-
-  output.close();
-
-  compileMEX("+" + basename, "static", mexext, { filename }, matlabroot, dynareroot);
-}
-
 void
 StaticModel::writeStaticJuliaFile(const string &basename) const
 {
@@ -985,7 +872,7 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool use_dll, c
       writeStaticBytecode(basename);
 
       if (use_dll)
-        writeStaticCFile(basename, mexext, matlabroot, dynareroot);
+        writeModelCFile<false>(basename, mexext, matlabroot, dynareroot);
       else if (julia)
         writeStaticJuliaFile(basename);
       else // M-files
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 5fa8bb974ad724257194fd9dd9913a40a61b108a..dcc222f00f67d2c0792966c37a3d35362a5b1f04 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -37,9 +37,6 @@ private:
   //! Writes static model file (standard Matlab version)
   void writeStaticMFile(const string &basename) const;
 
-  // Writes and compiles static model file (C version)
-  void writeStaticCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
-
   //! Writes static model file (Julia version)
   void writeStaticJuliaFile(const string &basename) const;