diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index cc011879ff4694c364ca526db88d54bc75f02c70..9e7723a5fee57e1e7237d17480fc9529335d360a 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -179,140 +179,6 @@ DynamicModel::additionalBlockTemporaryTerms(int blk,
     d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
 }
 
-int
-DynamicModel::nzeDeterministicJacobianForBlock(int blk) const
-{
-  BlockSimulationType simulation_type = blocks[blk].simulation_type;
-  int block_recursive_size = blocks[blk].getRecursiveSize();
-
-  int nze_deterministic = 0;
-  if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
-      || simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
-    nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(),
-                                 [=](const auto &kv) {
-                                   auto [eq, var, lag] = kv.first;
-                                   return eq >= block_recursive_size && var >= block_recursive_size;
-                                 });
-  else if (simulation_type == BlockSimulationType::solveBackwardSimple
-           || simulation_type == BlockSimulationType::solveForwardSimple
-           || simulation_type == BlockSimulationType::solveBackwardComplete
-           || simulation_type == BlockSimulationType::solveForwardComplete)
-    nze_deterministic = count_if(blocks_derivatives[blk].begin(), blocks_derivatives[blk].end(),
-                                 [=](const auto &kv) {
-                                   auto [eq, var, lag] = kv.first;
-                                   return lag == 0 && eq >= block_recursive_size && var >= block_recursive_size;
-                                 });
-  return nze_deterministic;
-}
-
-void
-DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
-{
-  temporary_terms_t temporary_terms; // Temp terms written so far
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    {
-      BlockSimulationType simulation_type = blocks[blk].simulation_type;
-      int block_size = blocks[blk].size;
-      int block_mfs_size = blocks[blk].mfs_size;
-
-      // Number of nonzero derivatives for the various Jacobians
-      int nze_stochastic = blocks_derivatives[blk].size();
-      int nze_deterministic = nzeDeterministicJacobianForBlock(blk);
-      int nze_other_endo = blocks_derivatives_other_endo[blk].size();
-      int nze_exo = blocks_derivatives_exo[blk].size();
-      int nze_exo_det = blocks_derivatives_exo_det[blk].size();
-
-      filesystem::path filename {packageDir(basename) / "+block" / ("dynamic_" + to_string(blk+1) + ".m")};
-      ofstream output{filename, ios::out | ios::binary};
-      if (!output.is_open())
-        {
-          cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-
-      output << "%" << endl
-             << "% " << filename.string() << " : Computes dynamic version of one block" << endl
-             << "%" << endl
-             << "% Warning : this file is generated automatically by Dynare" << endl
-             << "%           from model file (.mod)" << endl << endl
-             << "%" << endl;
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "function [y, T, g1, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode)" << endl;
-      else
-        output << "function [residual, y, T, g1, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode)" << endl;
-
-      output << "  % ////////////////////////////////////////////////////////////////////////" << endl
-             << "  % //" << "                     Block "sv.substr(static_cast<int>(log10(blk + 1))) << blk+1
-             << "                                        //" << endl
-             << "  % //                     Simulation type "
-             << BlockSim(simulation_type) << "  //" << endl
-             << "  % ////////////////////////////////////////////////////////////////////////" << endl;
-
-      if (simulation_type != BlockSimulationType::evaluateForward
-          && simulation_type != BlockSimulationType::evaluateBackward)
-        output << "  residual=zeros(" << block_mfs_size << ",1);" << endl;
-
-      output << "  if stochastic_mode" << endl
-             << "    g1_i=zeros(" << nze_stochastic << ",1);" << endl
-             << "    g1_j=zeros(" << nze_stochastic << ",1);" << endl
-             << "    g1_v=zeros(" << nze_stochastic << ",1);" << endl
-             << "    g1_x_i=zeros(" << nze_exo << ",1);" << endl
-             << "    g1_x_j=zeros(" << nze_exo << ",1);" << endl
-             << "    g1_x_v=zeros(" << nze_exo << ",1);" << endl
-             << "    g1_xd_i=zeros(" << nze_exo_det << ",1);" << endl
-             << "    g1_xd_j=zeros(" << nze_exo_det << ",1);" << endl
-             << "    g1_xd_v=zeros(" << nze_exo_det << ",1);" << endl
-             << "    g1_o_i=zeros(" << nze_other_endo << ",1);" << endl
-             << "    g1_o_j=zeros(" << nze_other_endo << ",1);" << endl
-             << "    g1_o_v=zeros(" << nze_other_endo << ",1);" << endl;
-      if (simulation_type != BlockSimulationType::evaluateForward
-          && simulation_type != BlockSimulationType::evaluateBackward)
-        output << "  else" << endl
-               << "    g1_i=zeros(" << nze_deterministic << ",1);" << endl
-               << "    g1_j=zeros(" << nze_deterministic << ",1);" << endl
-               << "    g1_v=zeros(" << nze_deterministic << ",1);" << endl;
-      output << "  end" << endl
-             << endl;
-
-      writeDynamicPerBlockHelper<ExprNodeOutputType::matlabDynamicModel>(blk, output, temporary_terms,
-                                 nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo);
-
-      output << endl
-             << "  if stochastic_mode" << endl
-             << "    g1=sparse(g1_i, g1_j, g1_v, " << block_size << ", " << blocks_jacob_cols_endo[blk].size() << ");" << endl
-             << "    varargout{1}=sparse(g1_x_i, g1_x_j, g1_x_v, " << block_size << ", " << blocks_jacob_cols_exo[blk].size() << ");" << endl
-             << "    varargout{2}=sparse(g1_xd_i, g1_xd_j, g1_xd_v, " << block_size << ", " << blocks_jacob_cols_exo_det[blk].size() << ");" << endl
-             << "    varargout{3}=sparse(g1_o_i, g1_o_j, g1_o_v, " << block_size << ", " << blocks_jacob_cols_other_endo[blk].size() << ");" << endl
-             << "  else" << endl;
-      switch (simulation_type)
-        {
-        case BlockSimulationType::evaluateForward:
-        case BlockSimulationType::evaluateBackward:
-          output << "    g1=[];" << endl;
-          break;
-        case BlockSimulationType::solveBackwardSimple:
-        case BlockSimulationType::solveForwardSimple:
-        case BlockSimulationType::solveBackwardComplete:
-        case BlockSimulationType::solveForwardComplete:
-          output << "    g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size
-                 << ", " << block_mfs_size << ");" << endl;
-          break;
-        case BlockSimulationType::solveTwoBoundariesSimple:
-        case BlockSimulationType::solveTwoBoundariesComplete:
-          output << "    g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size
-                 << ", " << 3*block_mfs_size << ");" << endl;
-          break;
-        default:
-          break;
-        }
-      output << "  end" << endl
-             << "end" << endl;
-      output.close();
-    }
-}
-
 void
 DynamicModel::writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block,
                                                       const temporary_terms_t &temporary_terms_union,
@@ -346,212 +212,6 @@ DynamicModel::writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file,
     }
 }
 
-vector<filesystem::path>
-DynamicModel::writeDynamicPerBlockCFiles(const string &basename, const string &mexext,
-                                         const filesystem::path &matlabroot,
-                                         const filesystem::path &dynareroot) const
-{
-  temporary_terms_t temporary_terms; // Temp terms written so far
-  const filesystem::path model_src_dir { filesystem::path{basename} / "model" / "src" };
-  vector<filesystem::path> compiled_object_files;
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    {
-      BlockSimulationType simulation_type = blocks[blk].simulation_type;
-      int block_size = blocks[blk].size;
-      int block_mfs_size = blocks[blk].mfs_size;
-
-      // Number of nonzero derivatives for the various Jacobians
-      int nze_stochastic = blocks_derivatives[blk].size();
-      int nze_deterministic = nzeDeterministicJacobianForBlock(blk);
-      int nze_other_endo = blocks_derivatives_other_endo[blk].size();
-      int nze_exo = blocks_derivatives_exo[blk].size();
-      int nze_exo_det = blocks_derivatives_exo_det[blk].size();
-
-      filesystem::path filename { model_src_dir / ("dynamic_" + to_string(blk+1) + ".c") };
-      ofstream output{filename, ios::out | ios::binary};
-      if (!output.is_open())
-        {
-          cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-
-      output << "/* Block " << blk+1 << endl
-             << "   " << BlockSim(simulation_type) << " */" << endl
-             << endl
-             << "#include <math.h>" << endl
-             << "#include <stdlib.h>" << endl
-             << "#include <stdbool.h>" << endl
-             << R"(#include "mex.h")" << endl
-             << endl;
-
-      // Write function definition if BinaryOpcode::powerDeriv is used
-      writePowerDerivHeader(output);
-
-      output << endl;
-
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "void dynamic_" << blk+1 << "(double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, double *restrict T, int it_, bool stochastic_mode, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v, double *restrict g1_x_i, double *restrict g1_x_j, double *restrict g1_x_v, double *restrict g1_xd_i, double *restrict g1_xd_j, double *restrict g1_xd_v, double *restrict g1_o_i, double *restrict g1_o_j, double *restrict g1_o_v)" << endl;
-      else
-        output << "void dynamic_" << blk+1 << "(double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, double *restrict T, int it_, bool stochastic_mode, double *restrict residual, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v, double *restrict g1_x_i, double *restrict g1_x_j, double *restrict g1_x_v, double *restrict g1_xd_i, double *restrict g1_xd_j, double *restrict g1_xd_v, double *restrict g1_o_i, double *restrict g1_o_j, double *restrict g1_o_v)" << endl;
-      output << '{' << endl;
-
-      writeDynamicPerBlockHelper<ExprNodeOutputType::CDynamicModel>(blk, output, temporary_terms,
-                                 nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo);
-
-      output << '}' << endl
-             << endl;
-
-      ostringstream header;
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        header << "void dynamic_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, const mxArray *steady_state, mxArray *T, const mxArray *it_, const mxArray *stochastic_mode, mxArray **g1, mxArray **g1_x, mxArray **g1_xd, mxArray **g1_o)";
-      else
-        header << "void dynamic_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, const mxArray *steady_state, mxArray *T, const mxArray *it_, const mxArray *stochastic_mode, mxArray **residual, mxArray **g1, mxArray **g1_x, mxArray **g1_xd, mxArray **g1_o)";
-      output << header.str() << endl
-             << '{' << endl
-             << "  int nb_row_x = mxGetM(x);" << endl;
-
-      if (simulation_type != BlockSimulationType::evaluateForward
-          && simulation_type != BlockSimulationType::evaluateBackward)
-        output << "  *residual = mxCreateDoubleMatrix(" << block_mfs_size << ",1,mxREAL);" << endl;
-
-      output << "  mxArray *g1_i = NULL, *g1_j = NULL, *g1_v = NULL;" << endl
-             << "  mxArray *g1_x_i = NULL, *g1_x_j = NULL, *g1_x_v = NULL;" << endl
-             << "  mxArray *g1_xd_i = NULL, *g1_xd_j = NULL, *g1_xd_v = NULL;" << endl
-             << "  mxArray *g1_o_i = NULL, *g1_o_j = NULL, *g1_o_v = NULL;" << endl
-             << "  if (mxGetScalar(stochastic_mode)) {" << endl
-             << "    g1_i=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
-             << "    g1_j=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
-             << "    g1_v=mxCreateDoubleMatrix(" << nze_stochastic << ",1,mxREAL);" << endl
-             << "    g1_x_i=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
-             << "    g1_x_j=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
-             << "    g1_x_v=mxCreateDoubleMatrix(" << nze_exo << ",1,mxREAL);" << endl
-             << "    g1_xd_i=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
-             << "    g1_xd_j=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
-             << "    g1_xd_v=mxCreateDoubleMatrix(" << nze_exo_det << ",1,mxREAL);" << endl
-             << "    g1_o_i=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl
-             << "    g1_o_j=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl
-             << "    g1_o_v=mxCreateDoubleMatrix(" << nze_other_endo << ",1,mxREAL);" << endl;
-      if (simulation_type != BlockSimulationType::evaluateForward
-          && simulation_type != BlockSimulationType::evaluateBackward)
-        output << "  } else {" << endl
-               << "    g1_i=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl
-               << "    g1_j=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl
-               << "    g1_v=mxCreateDoubleMatrix(" << nze_deterministic << ",1,mxREAL);" << endl;
-      output << "  }" << endl
-             << endl;
-
-      // N.B.: In the following, it_ is decreased by 1, to follow C convention
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "  dynamic_" << blk+1 << "(mxGetPr(y), mxGetPr(x), nb_row_x, mxGetPr(params), mxGetPr(steady_state), mxGetPr(T), mxGetScalar(it_)-1, mxGetScalar(stochastic_mode), g1_i ? mxGetPr(g1_i) : NULL, g1_j ? mxGetPr(g1_j) : NULL, g1_v ? mxGetPr(g1_v) : NULL, g1_x_i ? mxGetPr(g1_x_i) : NULL, g1_x_j ? mxGetPr(g1_x_j) : NULL, g1_x_v ? mxGetPr(g1_x_v) : NULL, g1_xd_i ? mxGetPr(g1_xd_i) : NULL, g1_xd_j ? mxGetPr(g1_xd_j) : NULL, g1_xd_v ? mxGetPr(g1_xd_v) : NULL, g1_o_i ? mxGetPr(g1_o_i) : NULL, g1_o_j ? mxGetPr(g1_o_j) : NULL, g1_o_v ? mxGetPr(g1_o_v) : NULL);" << endl;
-      else
-        output << "  dynamic_" << blk+1 << "(mxGetPr(y), mxGetPr(x), nb_row_x, mxGetPr(params), mxGetPr(steady_state), mxGetPr(T), mxGetScalar(it_)-1, mxGetScalar(stochastic_mode), mxGetPr(*residual), g1_i ? mxGetPr(g1_i) : NULL, g1_j ? mxGetPr(g1_j) : NULL, g1_v ? mxGetPr(g1_v) : NULL, g1_x_i ? mxGetPr(g1_x_i) : NULL, g1_x_j ? mxGetPr(g1_x_j) : NULL, g1_x_v ? mxGetPr(g1_x_v) : NULL, g1_xd_i ? mxGetPr(g1_xd_i) : NULL, g1_xd_j ? mxGetPr(g1_xd_j) : NULL, g1_xd_v ? mxGetPr(g1_xd_v) : NULL, g1_o_i ? mxGetPr(g1_o_i) : NULL, g1_o_j ? mxGetPr(g1_o_j) : NULL, g1_o_v ? mxGetPr(g1_o_v) : NULL);" << endl;
-
-      output << endl
-             << "  if (mxGetScalar(stochastic_mode)) {" << endl
-             << "    mxArray *m = mxCreateDoubleScalar(" << block_size << ");" << endl
-             << "    mxArray *n = mxCreateDoubleScalar(" << blocks_jacob_cols_endo[blk].size() << ");" << endl
-             << "    mxArray *plhs[1];" << endl
-             << "    mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
-             << R"(    mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
-             << "    *g1=plhs[0];" << endl
-             << "    mxDestroyArray(g1_i);" << endl
-             << "    mxDestroyArray(g1_j);" << endl
-             << "    mxDestroyArray(g1_v);" << endl
-             << "    mxDestroyArray(n);" << endl
-             << "    n = mxCreateDoubleScalar(" << blocks_jacob_cols_exo[blk].size() << ");" << endl
-             << "    mxArray *prhs_x[5] = { g1_x_i, g1_x_j, g1_x_v, m, n };" << endl
-             << R"(    mexCallMATLAB(1, plhs, 5, prhs_x, "sparse");)" << endl
-             << "    *g1_x=plhs[0];" << endl
-             << "    mxDestroyArray(g1_x_i);" << endl
-             << "    mxDestroyArray(g1_x_j);" << endl
-             << "    mxDestroyArray(g1_x_v);" << endl
-             << "    mxDestroyArray(n);" << endl
-             << "    n = mxCreateDoubleScalar(" << blocks_jacob_cols_exo_det[blk].size() << ");" << endl
-             << "    mxArray *prhs_xd[5] = { g1_xd_i, g1_xd_j, g1_xd_v, m, n };" << endl
-             << R"(    mexCallMATLAB(1, plhs, 5, prhs_xd, "sparse");)" << endl
-             << "    *g1_xd=plhs[0];" << endl
-             << "    mxDestroyArray(g1_xd_i);" << endl
-             << "    mxDestroyArray(g1_xd_j);" << endl
-             << "    mxDestroyArray(g1_xd_v);" << endl
-             << "    mxDestroyArray(n);" << endl
-             << "    n = mxCreateDoubleScalar(" << blocks_jacob_cols_other_endo[blk].size() << ");" << endl
-             << "    mxArray *prhs_o[5] = { g1_o_i, g1_o_j, g1_o_v, m, n };" << endl
-             << R"(    mexCallMATLAB(1, plhs, 5, prhs_o, "sparse");)" << endl
-             << "    *g1_o=plhs[0];" << endl
-             << "    mxDestroyArray(g1_o_i);" << endl
-             << "    mxDestroyArray(g1_o_j);" << endl
-             << "    mxDestroyArray(g1_o_v);" << endl
-             << "    mxDestroyArray(n);" << endl
-             << "    mxDestroyArray(m);" << endl
-             << "  } else {" << endl;
-      switch (simulation_type)
-        {
-        case BlockSimulationType::evaluateForward:
-        case BlockSimulationType::evaluateBackward:
-          output << "    *g1=mxCreateDoubleMatrix(0,0,mxREAL);" << endl;
-          break;
-        case BlockSimulationType::solveBackwardSimple:
-        case BlockSimulationType::solveForwardSimple:
-        case BlockSimulationType::solveBackwardComplete:
-        case BlockSimulationType::solveForwardComplete:
-          output << "    mxArray *m = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
-                 << "    mxArray *n = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
-                 << "    mxArray *plhs[1];" << endl
-                 << "    mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
-                 << R"(    mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
-                 << "    *g1=plhs[0];" << endl
-                 << "    mxDestroyArray(g1_i);" << endl
-                 << "    mxDestroyArray(g1_j);" << endl
-                 << "    mxDestroyArray(g1_v);" << endl
-                 << "    mxDestroyArray(n);" << endl
-                 << "    mxDestroyArray(m);" << endl;
-          break;
-        case BlockSimulationType::solveTwoBoundariesSimple:
-        case BlockSimulationType::solveTwoBoundariesComplete:
-          output << "    mxArray *m = mxCreateDoubleScalar(" << block_mfs_size << ");" << endl
-                 << "    mxArray *n = mxCreateDoubleScalar(" << 3*block_mfs_size << ");" << endl
-                 << "    mxArray *plhs[1];" << endl
-                 << "    mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
-                 << R"(    mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
-                 << "    *g1=plhs[0];" << endl
-                 << "    mxDestroyArray(g1_i);" << endl
-                 << "    mxDestroyArray(g1_j);" << endl
-                 << "    mxDestroyArray(g1_v);" << endl
-                 << "    mxDestroyArray(n);" << endl
-                 << "    mxDestroyArray(m);" << endl;
-          break;
-        default:
-          break;
-        }
-      output << "    *g1_x=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
-             << "    *g1_xd=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
-             << "    *g1_o=mxCreateDoubleMatrix(0,0,mxREAL);" << endl
-             << "  }" << endl
-             << "}" << endl;
-      output.close();
-
-      // Compile intermediary object under <MODFILE>/model/src/
-      compiled_object_files.emplace_back(compileMEX(model_src_dir, "dynamic_" + to_string(blk+1),
-                                                    mexext, { filename }, matlabroot, dynareroot,
-                                                    false));
-
-      filename = model_src_dir / ("dynamic_" + to_string(blk+1) + ".h");
-      ofstream header_output{filename, ios::out | ios::binary};
-      if (!header_output.is_open())
-        {
-          cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-      header_output << header.str() << ';' << endl;
-      header_output.close();
-    }
-  return compiled_object_files;
-}
-
 void
 DynamicModel::writeDynamicBytecode(const string &basename) const
 {
@@ -756,127 +416,6 @@ DynamicModel::printNonZeroHessianEquations(ostream &output) const
     output << "]";
 }
 
-void
-DynamicModel::writeDynamicBlockMFile(const string &basename) const
-{
-  filesystem::path filename {packageDir(basename) / "dynamic.m"};
-  ofstream output{filename, ios::out | ios::binary};
-  if (!output.is_open())
-    {
-      cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  output << "function [residual, y, T, g1, varargout] = dynamic(nblock, y, x, params, steady_state, T, it_, stochastic_mode)" << endl
-         << "  switch nblock" << endl;
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    {
-      output << "    case " << blk+1 << endl;
-
-      BlockSimulationType simulation_type = blocks[blk].simulation_type;
-
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "      [y, T, g1, varargout{1:nargout-4}] = " << basename << ".block.dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode);" << endl
-               << "      residual = [];" << endl;
-      else
-        output << "      [residual, y, T, g1, varargout{1:nargout-4}] = " << basename << ".block.dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, stochastic_mode);" << endl;
-    }
-  output << "  end" << endl
-         << "end" << endl;
-
-  output.close();
-}
-
-void
-DynamicModel::writeDynamicBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const
-{
-  const filesystem::path filename {basename + "/model/src/dynamic.c"};
-
-  ofstream output{filename, ios::out | ios::binary};
-  if (!output.is_open())
-    {
-      cerr << "Error: Can't open file " << filename.string() << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  output << "#include <math.h>" << endl
-         << R"(#include "mex.h")" << endl;
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    output << R"(#include "dynamic_)" << blk+1 << R"(.h")" << endl;
-
-  output << endl;
-  writePowerDeriv(output);
-
-  output << endl
-         << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
-         << "{" << endl
-         << "  if (nrhs != 8)" << endl
-         << R"(    mexErrMsgTxt("Requires exactly 8 input arguments");)" << endl
-         << "  if (nlhs > 7)" << endl
-         << R"(    mexErrMsgTxt("Accepts at most 7 output arguments");)" << endl
-         << "  int nblock = (int) mxGetScalar(prhs[0]);" << endl
-         << "  const mxArray *y = prhs[1], *x = prhs[2], *params = prhs[3], *steady_state = prhs[4], *T = prhs[5], *it_ = prhs[6], *stochastic_mode = prhs[7];" << endl
-         << "  mxArray *T_new = mxDuplicateArray(T);" << endl
-         << "  mxArray *y_new = mxDuplicateArray(y);" << endl
-         << "  mxArray *residual, *g1, *g1_x, *g1_xd, *g1_o;" << endl
-         << "  switch (nblock)" << endl
-         << "    {" << endl;
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    {
-      output << "    case " << blk+1 << ':' << endl;
-
-      BlockSimulationType simulation_type = blocks[blk].simulation_type;
-
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "      dynamic_" << blk+1 << "_mx(y_new, x, params, steady_state, T_new, it_, stochastic_mode, &g1, &g1_x, &g1_xd, &g1_o);" << endl
-               << "      residual = mxCreateDoubleMatrix(0,0,mxREAL);" << endl;
-      else
-        output << "      dynamic_" << blk+1 << "_mx(y_new, x, params, steady_state, T_new, it_, stochastic_mode, &residual, &g1, &g1_x, &g1_xd, &g1_o);" << endl;
-      output << "      break;" << endl;
-    }
-  output << "    }" << endl
-         << endl
-         << "  if (nlhs >= 1)" << endl
-         << "    plhs[0] = residual;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(residual);" << endl
-         << "  if (nlhs >= 2)" << endl
-         << "    plhs[1] = y_new;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(y_new);" << endl
-         << "  if (nlhs >= 3)" << endl
-         << "    plhs[2] = T_new;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(T_new);" << endl
-         << "  if (nlhs >= 4)" << endl
-         << "    plhs[3] = g1;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(g1);" << endl
-         << "  if (nlhs >= 5)" << endl
-         << "    plhs[4] = g1_x;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(g1_x);" << endl
-         << "  if (nlhs >= 6)" << endl
-         << "    plhs[5] = g1_xd;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(g1_xd);" << endl
-         << "  if (nlhs >= 7)" << endl
-         << "    plhs[6] = g1_o;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(g1_o);" << endl
-         << "}" << endl;
-
-  output.close();
-
-  per_block_object_files.push_back(filename);
-  compileMEX(packageDir(basename), "dynamic", mexext, per_block_object_files, matlabroot, dynareroot);
-}
-
 void
 DynamicModel::writeDynamicMWrapperFunction(const string &basename, const string &ending) const
 {
@@ -1353,8 +892,7 @@ DynamicModel::includeExcludeEquations(const string &inc_exc_option_value, bool e
 }
 
 void
-DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
-                                     const vector<int> &state_var, bool estimation_present) const
+DynamicModel::writeBlockDriverOutput(ostream &output) const
 {
   output << "M_.block_structure.time_recursive = " << boolalpha << time_recursive_block_decomposition << ";" << endl;
 
@@ -1362,14 +900,6 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
     {
       int block_size = blocks[blk].size;
       output << "M_.block_structure.block(" << blk+1 << ").Simulation_Type = " << static_cast<int>(blocks[blk].simulation_type) << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").maximum_lag = " << blocks[blk].max_lag << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").maximum_lead = " << blocks[blk].max_lead << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").maximum_endo_lag = " << blocks[blk].max_endo_lag << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").maximum_endo_lead = " << blocks[blk].max_endo_lead << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").maximum_exo_lag = " << blocks[blk].max_exo_lag << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").maximum_exo_lead = " << blocks[blk].max_exo_lead << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").maximum_exo_det_lag = " << blocks[blk].max_exo_det_lag << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").maximum_exo_det_lead = " << blocks[blk].max_exo_det_lead << ";" << endl
              << "M_.block_structure.block(" << blk+1 << ").endo_nbr = " << block_size << ";" << endl
              << "M_.block_structure.block(" << blk+1 << ").mfs = " << blocks[blk].mfs_size << ";" << endl
              << "M_.block_structure.block(" << blk+1 << ").equation = [";
@@ -1380,91 +910,6 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
       for (int var = 0; var < block_size; var++)
         output << " " << getBlockVariableID(blk, var)+1;
       output << "];" << endl
-             << "M_.block_structure.block(" << blk+1 << ").exogenous = [";
-      for (int exo : blocks_exo[blk])
-        output << " " << exo+1;
-      output << "];" << endl
-             << "M_.block_structure.block(" << blk+1 << ").exo_nbr = " << blocks_exo[blk].size() << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").exogenous_det = [";
-      for (int exo_det : blocks_exo_det[blk])
-        output << " " << exo_det+1;
-      output << "];" << endl
-             << "M_.block_structure.block(" << blk+1 << ").exo_det_nbr = " << blocks_exo_det[blk].size() << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").other_endogenous = [";
-      for (int other_endo : blocks_other_endo[blk])
-        output << " " << other_endo+1;
-      output << "];" << endl
-             << "M_.block_structure.block(" << blk+1 << ").other_endogenous_block = [";
-      for (int other_endo : blocks_other_endo[blk])
-        output << " " << endo2block[other_endo]+1;
-      output << "];" << endl;
-
-      output << "M_.block_structure.block(" << blk+1 << ").tm1 = zeros(" << blocks_other_endo[blk].size() << ", " << state_var.size() << ");" << endl;
-      for (int line{1};
-           auto other_endo : blocks_other_endo[blk])
-        {
-          if (auto it = find(state_var.begin(), state_var.end(), other_endo);
-              it != state_var.end())
-            output << "M_.block_structure.block(" << blk+1 << ").tm1("
-                   << line << ", "
-                   << distance(state_var.begin(), it)+1 << ") = 1;" << endl;
-          line++;
-        }
-
-      output << "M_.block_structure.block(" << blk+1 << ").other_endo_nbr = " << blocks_other_endo[blk].size() << ";" << endl;
-
-      int count_lead_lag_incidence = 0;
-      vector<int> local_state_var;
-      output << "M_.block_structure.block(" << blk+1 << ").lead_lag_incidence = [" << endl;
-      for (int lag = -1; lag <= 1; lag++)
-        {
-          for (int var = 0; var < block_size; var++)
-            {
-              for (int eq = 0; eq < block_size; eq++)
-                if (blocks_derivatives[blk].contains({ eq, var, lag }))
-                  {
-                    if (lag == -1)
-                      local_state_var.push_back(getBlockVariableID(blk, var));
-                    output << " " << ++count_lead_lag_incidence;
-                    goto var_found;
-                  }
-              output << " 0";
-            var_found:
-              ;
-            }
-          output << ";" << endl;
-        }
-      output << "];" << endl;
-
-      output << "M_.block_structure.block(" << blk+1 << ").sorted_col_dr_ghx = [";
-      for (int lsv : local_state_var)
-        output << distance(state_var.begin(), find(state_var.begin(), state_var.end(), lsv))+1 << " ";
-      output << "];" << endl;
-
-      count_lead_lag_incidence = 0;
-      output << "M_.block_structure.block(" << blk+1 << ").lead_lag_incidence_other = [" << endl;
-      for (int lag = -1; lag <= 1; lag++)
-        {
-          for (int other_endo : blocks_other_endo[blk])
-            {
-              for (int eq = 0; eq < block_size; eq++)
-                if (blocks_derivatives_other_endo[blk].contains({ eq, other_endo, lag }))
-                  {
-                    output << " " << ++count_lead_lag_incidence;
-                    goto other_endo_found;
-                  }
-              output << " 0";
-            other_endo_found:
-              ;
-            }
-          output << ";" << endl;
-        }
-      output << "];" << endl;
-
-      output << "M_.block_structure.block(" << blk+1 << ").n_static = " << blocks[blk].n_static << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").n_forward = " << blocks[blk].n_forward << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").n_backward = " << blocks[blk].n_backward << ";" << endl
-             << "M_.block_structure.block(" << blk+1 << ").n_mixed = " << blocks[blk].n_mixed << ";" << endl
              << "M_.block_structure.block(" << blk+1 << ").is_linear = " << boolalpha << blocks[blk].linear << ';' << endl
              << "M_.block_structure.block(" << blk+1 << ").NNZDerivatives = " << blocks_derivatives[blk].size() << ';' << endl;
     }
@@ -1499,138 +944,10 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
       output << "];" << endl;
     }
   output << "M_.block_structure.dyn_tmp_nbr = " << blocks_temporary_terms_idxs.size() << ';' << endl;
-
-  if (estimation_present)
-    {
-      filesystem::create_directories(basename + "/model/bytecode");
-      const filesystem::path main_name {basename + "/model/bytecode/kfi"};
-      ofstream KF_index_file{main_name, ios::out | ios::binary | ios::ate};
-      if (!KF_index_file.is_open())
-        {
-          cerr << "ERROR: Can't open file " << main_name.string() << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-      int n_obs = symbol_table.observedVariablesNbr();
-      int n_state = state_var.size();
-      for (int it : state_var)
-        if (symbol_table.isObservedVariable(symbol_table.getID(SymbolType::endogenous, it)))
-          n_obs--;
-
-      int n = n_obs + n_state;
-      output << "M_.nobs_non_statevar = " << n_obs << ";" << endl;
-      int nb_diag = 0;
-
-      vector<int> i_nz_state_var(n);
-      for (int i = 0; i < n_obs; i++)
-        i_nz_state_var[i] = n;
-      int lp = n_obs;
-
-      vector<int> state_equ;
-      for (int it : state_var)
-        state_equ.push_back(eq_idx_block2orig[endo_idx_orig2block[it]]);
-
-      for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-        {
-          int nze = 0;
-          for (int i = 0; i < blocks[blk].size; i++)
-            if (int var = getBlockVariableID(blk, i);
-                find(state_var.begin(), state_var.end(), var) != state_var.end())
-              nze++;
-
-          if (blk == 0)
-            {
-              set<pair<int, int>> row_state_var_incidence;
-              for (const auto &[idx, ignore] : blocks_derivatives[blk])
-                if (auto it_state_var = find(state_var.begin(), state_var.end(), getBlockVariableID(blk, get<1>(idx)));
-                    it_state_var != state_var.end())
-                  if (auto it_state_equ = find(state_equ.begin(), state_equ.end(), getBlockEquationID(blk, get<0>(idx)));
-                      it_state_equ != state_equ.end())
-                    row_state_var_incidence.emplace(it_state_equ - state_equ.begin(), it_state_var - state_var.begin());
-              auto row_state_var_incidence_it = row_state_var_incidence.begin();
-              bool diag = true;
-              int nb_diag_r = 0;
-              while (row_state_var_incidence_it != row_state_var_incidence.end() && diag)
-                {
-                  diag = (row_state_var_incidence_it->first == row_state_var_incidence_it->second);
-                  if (diag)
-                    {
-                      int equ = row_state_var_incidence_it->first;
-                      row_state_var_incidence_it++;
-                      if (equ != row_state_var_incidence_it->first)
-                        nb_diag_r++;
-                    }
-
-                }
-              set<pair<int, int>> col_state_var_incidence;
-              for (auto [equ, var] : row_state_var_incidence)
-                col_state_var_incidence.emplace(var, equ);
-              auto col_state_var_incidence_it = col_state_var_incidence.begin();
-              diag = true;
-              int nb_diag_c = 0;
-              while (col_state_var_incidence_it != col_state_var_incidence.end() && diag)
-                {
-                  diag = (col_state_var_incidence_it->first == col_state_var_incidence_it->second);
-                  if (diag)
-                    {
-                      int var = col_state_var_incidence_it->first;
-                      col_state_var_incidence_it++;
-                      if (var != col_state_var_incidence_it->first)
-                        nb_diag_c++;
-                    }
-                }
-              nb_diag = min(nb_diag_r, nb_diag_c);
-              row_state_var_incidence.clear();
-              col_state_var_incidence.clear();
-            }
-          for (int i = 0; i < nze; i++)
-            i_nz_state_var[lp + i] = lp + nze;
-          lp += nze;
-        }
-      output << "M_.nz_state_var = [";
-      for (int i = 0; i < lp; i++)
-        output << i_nz_state_var[i] << " ";
-      output << "];" << endl
-             << "M_.n_diag = " << nb_diag << ";" << endl;
-      KF_index_file.write(reinterpret_cast<char *>(&nb_diag), sizeof(nb_diag));
-
-      using index_KF = pair<int, pair<int, int >>;
-      vector<index_KF> v_index_KF;
-      for (int i = 0; i < n; i++)
-        for (int j = n_obs; j < n; j++)
-          {
-            int j1 = j - n_obs;
-            int j1_n_state = j1 * n_state - n_obs;
-            if ((i < n_obs) || (i >= nb_diag + n_obs) || (j1 >= nb_diag))
-              for (int k = n_obs; k < i_nz_state_var[i]; k++)
-                v_index_KF.emplace_back(i + j1 * n, pair(i + k * n, k + j1_n_state));
-          }
-      int size_v_index_KF = v_index_KF.size();
-
-      KF_index_file.write(reinterpret_cast<char *>(&size_v_index_KF), sizeof(size_v_index_KF));
-      for (auto &it : v_index_KF)
-        KF_index_file.write(reinterpret_cast<char *>(&it), sizeof(index_KF));
-
-      vector<index_KF> v_index_KF_2;
-      int n_n_obs = n * n_obs;
-      for (int i = 0; i < n; i++)
-        for (int j = i; j < n; j++)
-          if ((i < n_obs) || (i >= nb_diag + n_obs) || (j < n_obs) || (j >= nb_diag + n_obs))
-            for (int k = n_obs; k < i_nz_state_var[j]; k++)
-              {
-                int k_n = k * n;
-                v_index_KF_2.emplace_back(i * n + j, pair(i + k_n - n_n_obs, j + k_n));
-              }
-      int size_v_index_KF_2 = v_index_KF_2.size();
-
-      KF_index_file.write(reinterpret_cast<char *>(&size_v_index_KF_2), sizeof(size_v_index_KF_2));
-      for (auto &it : v_index_KF_2)
-        KF_index_file.write(reinterpret_cast<char *>(&it), sizeof(index_KF));
-      KF_index_file.close();
-    }
 }
 
 void
-DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool estimation_present, bool compute_xrefs) const
+DynamicModel::writeDriverOutput(ostream &output, bool compute_xrefs) const
 {
   /* Writing initialisation for M_.lead_lag_incidence matrix
      M_.lead_lag_incidence is a matrix with as many columns as there are
@@ -1760,7 +1077,7 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool es
 
   // Write the block structure of the model
   if (block_decomposed)
-    writeBlockDriverOutput(output, basename, state_var, estimation_present);
+    writeBlockDriverOutput(output);
 
   output << "M_.state_var = [";
   for (int it : state_var)
@@ -3321,32 +2638,12 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool use_dll,
   create_directories(model_dir / "bytecode" / "block");
 
   // Legacy representation
-  if (block)
-    {
-      if (use_dll)
-        {
-          auto per_block_object_files { writeDynamicPerBlockCFiles(basename, mexext, matlabroot, dynareroot) };
-          writeDynamicBlockCFile(basename, move(per_block_object_files), mexext, matlabroot, dynareroot);
-        }
-      else if (julia)
-        {
-          cerr << "'block' option is not available with Julia" << endl;
-          exit(EXIT_FAILURE);
-        }
-      else // M-files
-        {
-          writeDynamicPerBlockMFiles(basename);
-          writeDynamicBlockMFile(basename);
-        }
-    }
-  else
-    {
-      if (use_dll)
-        writeModelCFile<true>(basename, mexext, matlabroot, dynareroot);
-      else if (!julia) // M-files
-        writeDynamicMFile(basename);
-      // The legacy representation is no longer produced for Julia
-    }
+  if (use_dll)
+    writeModelCFile<true>(basename, mexext, matlabroot, dynareroot);
+  else if (!julia) // M-files
+    writeDynamicMFile(basename);
+  // The legacy representation is no longer produced for Julia
+
   writeDynamicBytecode(basename);
   if (block_decomposed)
     writeDynamicBlockBytecode(basename);
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index a71b08cc222cc9e36db2ee1c19d03d48c951de37..d8f04dcb3140bd761fa1688f5714f7496946bfcf 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -120,26 +120,6 @@ private:
 
   // Writes dynamic model file (MATLAB/Octave version, legacy representation)
   void writeDynamicMFile(const string &basename) const;
-  /* Writes the main dynamic function of block decomposed model (MATLAB/Octave
-     version, legacy representation) */
-  void writeDynamicBlockMFile(const string &basename) const;
-  /* Writes the main dynamic functions of block decomposed model (C version,
-     legacy representation), then compiles it with the per-block functions into
-     a single MEX */
-  void writeDynamicBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
-  /* Computes the number of nonzero elements in deterministic Jacobian of
-     block-decomposed model */
-  int nzeDeterministicJacobianForBlock(int blk) const;
-  // Helper for writing the per-block dynamic files of block decomposed models (legacy representation)
-  template<ExprNodeOutputType output_type>
-  void writeDynamicPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const;
-  /* Writes the per-block dynamic files of block decomposed model (MATLAB/Octave
-     version, legacy representation) */
-  void writeDynamicPerBlockMFiles(const string &basename) const;
-  /* Writes the per-block dynamic files of block decomposed model (C version,
-     legacy representation).
-     Returns the list of paths to the generated C source files (not the headers) */
-  vector<filesystem::path> writeDynamicPerBlockCFiles(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
   //! Writes the code of the block-decomposed model in virtual machine bytecode
   void writeDynamicBlockBytecode(const string &basename) const;
   // Writes derivatives w.r.t. exo, exo det and other endogenous
@@ -153,8 +133,7 @@ private:
   void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const;
 
   // Write the block structure of the model in the driver file
-  void writeBlockDriverOutput(ostream &output, const string &basename,
-                              const vector<int> &state_var, bool estimation_present) const;
+  void writeBlockDriverOutput(ostream &output) const;
 
   // Used by determineBlockDerivativesType()
   enum class BlockDerivativeType
@@ -334,7 +313,7 @@ public:
   void computingPass(int derivsOrder, int paramsDerivsOrder, const eval_context_t &eval_context,
                      bool no_tmp_terms, bool block, bool use_dll);
   //! Writes information about the dynamic model to the driver file
-  void writeDriverOutput(ostream &output, const string &basename, bool estimation_present, bool compute_xrefs) const;
+  void writeDriverOutput(ostream &output, bool compute_xrefs) const;
 
   //! Write JSON AST
   void writeJsonAST(ostream &output) const;
@@ -697,171 +676,6 @@ public:
   set<int> findPacExpectationEquationNumbers() const;
 };
 
-template<ExprNodeOutputType output_type>
-void
-DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms,
-                                         int nze_stochastic, int nze_deterministic, int nze_exo,
-                                         int nze_exo_det, int nze_other_endo) const
-{
-  static_assert(!isSparseModelOutput(output_type));
-
-  BlockSimulationType simulation_type { blocks[blk].simulation_type };
-  int block_mfs_size { blocks[blk].mfs_size };
-  int block_recursive_size { blocks[blk].getRecursiveSize() };
-
-  // Write residuals and temporary terms (incl. for derivatives)
-  writePerBlockHelper<output_type>(blk, output, temporary_terms);
-
-  if constexpr(isCOutput(output_type))
-    output << "  if (stochastic_mode) {" << endl;
-  else
-    output << "  if stochastic_mode" << endl;
-
-  ostringstream i_output, j_output, v_output;
-  int line_counter { ARRAY_SUBSCRIPT_OFFSET(output_type) };
-  for (const auto &[indices, d] : blocks_derivatives[blk])
-    {
-      const auto &[eq, var, lag] {indices};
-      int jacob_col { blocks_jacob_cols_endo[blk].at({ var, lag }) };
-      i_output << "    g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
-      j_output << "    g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
-      v_output << "    g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
-      d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
-      v_output << ';' << endl;
-      line_counter++;
-    }
-  assert(line_counter == nze_stochastic+ARRAY_SUBSCRIPT_OFFSET(output_type));
-  output << i_output.str() << j_output.str() << v_output.str();
-
-  i_output.str("");
-  j_output.str("");
-  v_output.str("");
-  line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
-  for (const auto &[indices, d] : blocks_derivatives_exo[blk])
-    {
-      const auto &[eq, var, lag] {indices};
-      int jacob_col { blocks_jacob_cols_exo[blk].at({ var, lag }) };
-      i_output << "    g1_x_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
-      j_output << "    g1_x_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
-      v_output << "    g1_x_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
-      d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
-      v_output << ';' << endl;
-      line_counter++;
-    }
-  assert(line_counter == nze_exo+ARRAY_SUBSCRIPT_OFFSET(output_type));
-  output << i_output.str() << j_output.str() << v_output.str();
-
-  i_output.str("");
-  j_output.str("");
-  v_output.str("");
-  line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
-  for (const auto &[indices, d] : blocks_derivatives_exo_det[blk])
-    {
-      const auto &[eq, var, lag] {indices};
-      int jacob_col { blocks_jacob_cols_exo_det[blk].at({ var, lag }) };
-      i_output << "    g1_xd_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
-      j_output << "    g1_xd_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
-      v_output << "    g1_xd_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
-      d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
-      v_output << ';' << endl;
-      line_counter++;
-    }
-  assert(line_counter == nze_exo_det+ARRAY_SUBSCRIPT_OFFSET(output_type));
-  output << i_output.str() << j_output.str() << v_output.str();
-
-  i_output.str("");
-  j_output.str("");
-  v_output.str("");
-  line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
-  for (const auto &[indices, d] : blocks_derivatives_other_endo[blk])
-    {
-      const auto &[eq, var, lag] {indices};
-      int jacob_col { blocks_jacob_cols_other_endo[blk].at({ var, lag }) };
-      i_output << "    g1_o_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
-      j_output << "    g1_o_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
-      v_output << "    g1_o_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-               << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
-      d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
-      v_output << ';' << endl;
-      line_counter++;
-    }
-  assert(line_counter == nze_other_endo+ARRAY_SUBSCRIPT_OFFSET(output_type));
-  output << i_output.str() << j_output.str() << v_output.str();
-
-  // Deterministic mode
-  if (simulation_type != BlockSimulationType::evaluateForward
-      && simulation_type != BlockSimulationType::evaluateBackward)
-    {
-      if constexpr(isCOutput(output_type))
-        output << "  } else {" << endl;
-      else
-        output << "  else" << endl;
-      i_output.str("");
-      j_output.str("");
-      v_output.str("");
-      line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
-      if (simulation_type == BlockSimulationType::solveBackwardSimple
-          || simulation_type == BlockSimulationType::solveForwardSimple
-          || simulation_type == BlockSimulationType::solveBackwardComplete
-          || simulation_type == BlockSimulationType::solveForwardComplete)
-        for (const auto &[indices, d] : blocks_derivatives[blk])
-          {
-            const auto &[eq, var, lag] {indices};
-            if (lag == 0 && eq >= block_recursive_size && var >= block_recursive_size)
-              {
-                i_output << "    g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                         << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
-                         << eq+1-block_recursive_size << ';' << endl;
-                j_output << "    g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                         << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
-                         << var+1-block_recursive_size << ';' << endl;
-                v_output << "    g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                         << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
-                d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
-                v_output << ';' << endl;
-                line_counter++;
-              }
-          }
-      else // solveTwoBoundariesSimple || solveTwoBoundariesComplete
-        for (const auto &[indices, d] : blocks_derivatives[blk])
-        {
-          const auto &[eq, var, lag] {indices};
-          assert(lag >= -1 && lag <= 1);
-          if (eq >= block_recursive_size && var >= block_recursive_size)
-            {
-              i_output << "    g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
-                       << eq+1-block_recursive_size << ';' << endl;
-              j_output << "    g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
-                       << var+1-block_recursive_size+block_mfs_size*(lag+1) << ';' << endl;
-              v_output << "    g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
-              d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
-              v_output << ';' << endl;
-              line_counter++;
-            }
-        }
-      assert(line_counter == nze_deterministic+ARRAY_SUBSCRIPT_OFFSET(output_type));
-      output << i_output.str() << j_output.str() << v_output.str();
-    }
-  if constexpr(isCOutput(output_type))
-    output << "  }" << endl;
-  else
-    output << "  end" << endl;
-}
-
 template<bool julia>
 void
 DynamicModel::writeParamsDerivativesFile(const string &basename) const
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 4e7c175f561fabd20ff98a6240537efda8e6d426..39212e7fecb0c68c04bb4215547cc90ea6b61dad 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -367,16 +367,6 @@ ModFile::checkPass(bool nostrict, bool stochastic)
           exit(EXIT_FAILURE);
         }
     }
-
-  // See dynare#1726
-  if ((stochastic_statement_present || mod_file_struct.check_present) && block && dynamic_model.mfs > 0)
-    {
-      /* NB: If mfs>0 but “block” is not passed, the block-DR routines will not
-         be called, so do not fail in that case (we may want to use the sparse
-         block representation) */
-      cerr << "ERROR: the `block` option used in conjunction with `mfs > 0` is incompatible with check, stoch_simul, estimation, osr, ramsey_policy, discretionary_policy, calib_smoother, identification, methods_of_moments and sensitivity commands" << endl;
-      exit(EXIT_FAILURE);
-    }
 }
 
 void
@@ -946,7 +936,7 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
 
   if (dynamic_model.equation_number() > 0)
     {
-      dynamic_model.writeDriverOutput(mOutputFile, basename, mod_file_struct.estimation_present, compute_xrefs);
+      dynamic_model.writeDriverOutput(mOutputFile, compute_xrefs);
       if (!no_static)
         static_model.writeDriverOutput(mOutputFile);
     }
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 01e401d5a5f8a1cf5ec852d2852285eb380e386d..2ac05d7ee626bcb7151445b9f879836c34cc7212 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -327,7 +327,6 @@ protected:
   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)
-  // This part is common to the legacy and sparse representations
   template<ExprNodeOutputType output_type>
   void writePerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
 
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index faa81502678be041a619955727f544f337d86470..24f7f5fdb7572417f9e1db6126c4b9d4fdf7d952 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -95,160 +95,6 @@ StaticModel::StaticModel(const DynamicModel &m) :
   user_set_compiler = m.user_set_compiler;
 }
 
-void
-StaticModel::writeStaticPerBlockMFiles(const string &basename) const
-{
-  temporary_terms_t temporary_terms; // Temp terms written so far
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    {
-      BlockSimulationType simulation_type = blocks[blk].simulation_type;
-
-      filesystem::path filename {packageDir(basename) / "+block" / ("static_" + to_string(blk+1) + ".m")};
-      ofstream output{filename, ios::out | ios::binary};
-      if (!output.is_open())
-        {
-          cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-      output << "%" << endl
-             << "% " << filename.string() << " : Computes static version of one block" << endl
-             << "%" << endl
-             << "% Warning : this file is generated automatically by Dynare" << endl
-             << "%           from model file (.mod)" << endl << endl
-             << "%" << endl;
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "function [y, T] = static_" << blk+1 << "(y, x, params, T)" << endl;
-      else
-        output << "function [residual, y, T, g1] = static_" << blk+1 << "(y, x, params, T)" << endl;
-
-      output << "  % ////////////////////////////////////////////////////////////////////////" << endl
-             << "  % //" << "                     Block "sv.substr(static_cast<int>(log10(blk + 1))) << blk+1
-             << "                                        //" << endl
-             << "  % //                     Simulation type "
-             << BlockSim(simulation_type) << "  //" << endl
-             << "  % ////////////////////////////////////////////////////////////////////////" << endl;
-
-      if (simulation_type != BlockSimulationType::evaluateBackward
-          && simulation_type != BlockSimulationType::evaluateForward)
-        output << "  residual=zeros(" << blocks[blk].mfs_size << ",1);" << endl
-               << "  g1_i=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
-               << "  g1_j=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
-               << "  g1_v=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
-               << endl;
-
-      writeStaticPerBlockHelper<ExprNodeOutputType::matlabStaticModel>(blk, output, temporary_terms);
-
-      if (simulation_type != BlockSimulationType::evaluateBackward
-          && simulation_type != BlockSimulationType::evaluateForward)
-        output << endl
-               << "  g1=sparse(g1_i, g1_j, g1_v, "  << blocks[blk].mfs_size << "," << blocks[blk].mfs_size << ");" << endl;
-
-      output << "end" << endl;
-      output.close();
-    }
-}
-
-vector<filesystem::path>
-StaticModel::writeStaticPerBlockCFiles(const string &basename, const string &mexext,
-                                       const filesystem::path &matlabroot,
-                                       const filesystem::path &dynareroot) const
-{
-  temporary_terms_t temporary_terms; // Temp terms written so far
-  const filesystem::path model_src_dir { filesystem::path{basename} / "model" / "src" };
-  vector<filesystem::path> compiled_object_files;
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    {
-      BlockSimulationType simulation_type = blocks[blk].simulation_type;
-
-      filesystem::path filename { model_src_dir / ("static_" + to_string(blk+1) + ".c") };
-      ofstream output{filename, ios::out | ios::binary};
-      if (!output.is_open())
-        {
-          cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-      output << "/* Block " << blk+1 << endl
-             << "   " << BlockSim(simulation_type) << " */" << endl
-             << endl
-             << "#include <math.h>" << endl
-             << "#include <stdlib.h>" << endl
-             << R"(#include "mex.h")" << endl
-             << endl;
-
-      // Write function definition if BinaryOpcode::powerDeriv is used
-      writePowerDerivHeader(output);
-
-      output << endl;
-
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "void static_" << blk+1 << "(double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl;
-      else
-        output << "void static_" << blk+1 << "(double *restrict y, const double *restrict x, const double *restrict params, double *restrict T, double *restrict residual, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v)" << endl;
-      output << '{' << endl;
-
-      writeStaticPerBlockHelper<ExprNodeOutputType::CStaticModel>(blk, output, temporary_terms);
-
-      output << '}' << endl
-             << endl;
-
-      ostringstream header;
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        {
-          header << "void static_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, mxArray *T)";
-          output << header.str() << endl
-                 << '{' << endl
-                 << "  static_" << blk+1 << "(mxGetPr(y), mxGetPr(x), mxGetPr(params), mxGetPr(T));" << endl
-                 << '}' << endl;
-        }
-      else
-        {
-          header << "void static_" << blk+1 << "_mx(mxArray *y, const mxArray *x, const mxArray *params, mxArray *T, mxArray **residual, mxArray **g1)";
-          output << header.str() << endl
-                 << '{' << endl
-                 << "  *residual = mxCreateDoubleMatrix(" << blocks[blk].mfs_size << ",1,mxREAL);" << endl
-                 << "  mxArray *g1_i = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
-                 << "  mxArray *g1_j = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
-                 << "  mxArray *g1_v = mxCreateDoubleMatrix(" << blocks_derivatives[blk].size() << ",1,mxREAL);" << endl
-                 << "  static_" << blk+1 << "(mxGetPr(y), mxGetPr(x), mxGetPr(params), mxGetPr(T), mxGetPr(*residual), mxGetPr(g1_i), mxGetPr(g1_j), mxGetPr(g1_v));" << endl
-                 << "  mxArray *plhs[1];" << endl
-                 << "  mxArray *m = mxCreateDoubleScalar(" << blocks[blk].mfs_size << ");" << endl
-                 << "  mxArray *n = mxCreateDoubleScalar(" << blocks[blk].mfs_size << ");" << endl
-                 << "  mxArray *prhs[5] = { g1_i, g1_j, g1_v, m, n };" << endl
-                 << R"(  mexCallMATLAB(1, plhs, 5, prhs, "sparse");)" << endl
-                 << "  *g1 = plhs[0];" << endl
-                 << "  mxDestroyArray(g1_i);" << endl
-                 << "  mxDestroyArray(g1_j);" << endl
-                 << "  mxDestroyArray(g1_v);" << endl
-                 << "  mxDestroyArray(m);" << endl
-                 << "  mxDestroyArray(n);" << endl
-                 << '}' << endl;
-        }
-
-      output.close();
-
-      // Compile intermediary object under <MODFILE>/model/src/
-      compiled_object_files.emplace_back(compileMEX(model_src_dir, "static_" + to_string(blk+1),
-                                                    mexext, { filename }, matlabroot, dynareroot,
-                                                    false));
-
-      filename = model_src_dir / ("static_" + to_string(blk+1) + ".h");
-      ofstream header_output{filename, ios::out | ios::binary};
-      if (!header_output.is_open())
-        {
-          cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-          exit(EXIT_FAILURE);
-        }
-      header_output << header.str() << ';' << endl;
-      header_output.close();
-    }
-  return compiled_object_files;
-}
-
 void
 StaticModel::writeStaticBytecode(const string &basename) const
 {
@@ -645,32 +491,12 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool use_dll, c
   create_directories(model_dir / "bytecode" / "block");
 
   // Legacy representation
-  if (block)
-    {
-      if (use_dll)
-        {
-          auto per_block_object_files { writeStaticPerBlockCFiles(basename, mexext, matlabroot, dynareroot) };
-          writeStaticBlockCFile(basename, move(per_block_object_files), mexext, matlabroot, dynareroot);
-        }
-      else if (julia)
-        {
-          cerr << "'block' option is not available with Julia" << endl;
-          exit(EXIT_FAILURE);
-        }
-      else // M-files
-        {
-          writeStaticPerBlockMFiles(basename);
-          writeStaticBlockMFile(basename);
-        }
-    }
-  else
-    {
-      if (use_dll)
-        writeModelCFile<false>(basename, mexext, matlabroot, dynareroot);
-      else if (!julia) // M-files
-        writeStaticMFile(basename);
-      // The legacy representation is no longer produced for Julia
-    }
+  if (use_dll)
+    writeModelCFile<false>(basename, mexext, matlabroot, dynareroot);
+  else if (!julia) // M-files
+    writeStaticMFile(basename);
+  // The legacy representation is no longer produced for Julia
+
   writeStaticBytecode(basename);
   if (block_decomposed)
     writeStaticBlockBytecode(basename);
@@ -699,117 +525,6 @@ StaticModel::exoPresentInEqs() const
   return false;
 }
 
-void
-StaticModel::writeStaticBlockMFile(const string &basename) const
-{
-  filesystem::path filename {packageDir(basename) / "static.m"};
-
-  ofstream output{filename, ios::out | ios::binary};
-  if (!output.is_open())
-    {
-      cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  output << "function [residual, y, T, g1] = static(nblock, y, x, params, T)" << endl
-         << "  switch nblock" << endl;
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    {
-      output << "    case " << blk+1 << endl;
-
-      BlockSimulationType simulation_type = blocks[blk].simulation_type;
-
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "      [y, T] = " << basename << ".block.static_" << blk+1 << "(y, x, params, T);" << endl
-               << "      residual = [];" << endl
-               << "      g1 = [];" << endl;
-      else
-        output << "      [residual, y, T, g1] = " << basename << ".block.static_" << blk+1 << "(y, x, params, T);" << endl;
-
-    }
-  output << "  end" << endl
-         << "end" << endl;
-  output.close();
-}
-
-void
-StaticModel::writeStaticBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const
-{
-  const filesystem::path filename {basename + "/model/src/static.c"};
-
-  ofstream output{filename, ios::out | ios::binary};
-  if (!output.is_open())
-    {
-      cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
-      exit(EXIT_FAILURE);
-    }
-
-  output << "#include <math.h>" << endl
-         << R"(#include "mex.h")" << endl;
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    output << R"(#include "static_)" << blk+1 << R"(.h")" << endl;
-
-  output << endl;
-  writePowerDeriv(output);
-
-  output << endl
-         << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
-         << "{" << endl
-         << "  if (nrhs != 5)" << endl
-         << R"(    mexErrMsgTxt("Requires exactly 5 input arguments");)" << endl
-         << "  if (nlhs > 4)" << endl
-         << R"(    mexErrMsgTxt("Accepts at most 4 output arguments");)" << endl
-         << "  int nblock = (int) mxGetScalar(prhs[0]);" << endl
-         << "  const mxArray *y = prhs[1], *x = prhs[2], *params = prhs[3], *T = prhs[4];" << endl
-         << "  mxArray *T_new = mxDuplicateArray(T);" << endl
-         << "  mxArray *y_new = mxDuplicateArray(y);" << endl
-         << "  mxArray *residual, *g1;" << endl
-         << "  switch (nblock)" << endl
-         << "    {" << endl;
-
-  for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
-    {
-      output << "    case " << blk+1 << ':' << endl;
-
-      BlockSimulationType simulation_type = blocks[blk].simulation_type;
-
-      if (simulation_type == BlockSimulationType::evaluateBackward
-          || simulation_type == BlockSimulationType::evaluateForward)
-        output << "      static_" << blk+1 << "_mx(y_new, x, params, T_new);" << endl
-               << "      residual = mxCreateDoubleMatrix(0,0,mxREAL);" << endl
-               << "      g1 = mxCreateDoubleMatrix(0,0,mxREAL);" << endl;
-      else
-        output << "      static_" << blk+1 << "_mx(y_new, x, params, T_new, &residual, &g1);" << endl;
-      output << "      break;" << endl;
-    }
-  output << "    }" << endl
-         << endl
-         << "  if (nlhs >= 1)" << endl
-         << "    plhs[0] = residual;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(residual);" << endl
-         << "  if (nlhs >= 2)" << endl
-         << "    plhs[1] = y_new;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(y_new);" << endl
-         << "  if (nlhs >= 3)" << endl
-         << "    plhs[2] = T_new;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(T_new);" << endl
-         << "  if (nlhs >= 4)" << endl
-         << "    plhs[3] = g1;" << endl
-         << "  else" << endl
-         << "    mxDestroyArray(g1);" << endl
-         << "}" << endl;
-  output.close();
-
-  per_block_object_files.push_back(filename);
-  compileMEX(packageDir(basename), "static", mexext, per_block_object_files, matlabroot, dynareroot);
-}
-
 void
 StaticModel::writeDriverOutput(ostream &output) const
 {
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 3d4ebd06fde41e33832956c727d3a560e4ed082a..50c73b8ec54380353fb98a785bcb1c733b4eb777 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2003-2022 Dynare Team
+ * Copyright © 2003-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -37,28 +37,6 @@ private:
   // Writes static model file (MATLAB/Octave version, legacy representation)
   void writeStaticMFile(const string &basename) const;
 
-  /* Writes the main static function of block decomposed model (MATLAB/Octave
-     version, legacy representation) */
-  void writeStaticBlockMFile(const string &basename) const;
-
-  /* Writes the main static functions of block decomposed model (C version,
-     legacy representation), then compiles it with the per-block functions into
-     a single MEX */
-  void writeStaticBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
-
-  // Helper for writing a per-block static file of block decomposed model (legacy representation)
-  template<ExprNodeOutputType output_type>
-  void writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
-
-  /* Writes the per-block static files of block decomposed model (MATLAB/Octave
-     version, legacy representation) */
-  void writeStaticPerBlockMFiles(const string &basename) const;
-
-  /* Writes the per-block static files of block decomposed model (C version,
-     legacy representation).
-     Returns the list of paths to the generated C source files (not the headers) */
-  vector<filesystem::path> writeStaticPerBlockCFiles(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
-
   //! Writes the code of the block-decomposed model in virtual machine bytecode
   void writeStaticBlockBytecode(const string &basename) const;
 
@@ -181,43 +159,6 @@ public:
   void addAllParamDerivId(set<int> &deriv_id_set) override;
 };
 
-template<ExprNodeOutputType output_type>
-void
-StaticModel::writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
-{
-  static_assert(!isSparseModelOutput(output_type));
-
-  BlockSimulationType simulation_type { blocks[blk].simulation_type };
-  int block_recursive_size { blocks[blk].getRecursiveSize() };
-
-  // Write residuals and temporary terms (incl. for derivatives)
-  writePerBlockHelper<output_type>(blk, output, temporary_terms);
-
-  // The Jacobian if we have to solve the block
-  if (simulation_type != BlockSimulationType::evaluateBackward
-      && simulation_type != BlockSimulationType::evaluateForward)
-    {
-      ostringstream i_output, j_output, v_output;
-      for (int line_counter { ARRAY_SUBSCRIPT_OFFSET(output_type) };
-           const auto &[indices, d] : blocks_derivatives[blk])
-        {
-          const auto &[eq, var, ignore] {indices};
-          i_output << "  g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1-block_recursive_size
-                   << ';' << endl;
-          j_output << "  g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << var+1-block_recursive_size
-                   << ';' << endl;
-          v_output << "  g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
-                   << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
-          d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
-          v_output << ';' << endl;
-          line_counter++;
-        }
-      output << i_output.str() << j_output.str() << v_output.str();
-    }
-}
-
 template<bool julia>
 void
 StaticModel::writeParamsDerivativesFile(const string &basename) const