diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 305d0de0..1f29299d 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -227,7 +227,7 @@ void
 DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
 {
   temporary_terms_t temporary_terms; // Temp terms written so far
-  constexpr ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabDynamicBlockModel;
+  constexpr ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabDynamicModel;
 
   for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
     {
@@ -1294,10 +1294,10 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
                     << "function [residual, dr] = dynamic(options_, M_, oo_, y, x, params, steady_state, it_, dr)" << endl;
 
   if (blocks_temporary_terms_idxs.size() > 0)
-    mDynamicModelFile << "  T=NaN(" << blocks_temporary_terms_idxs.size()
-                      << ",options_.periods+M_.maximum_lag+M_.maximum_lead);" << endl;
+    mDynamicModelFile << "  T=NaN(" << blocks_temporary_terms_idxs.size() << ");" << endl;
 
-  mDynamicModelFile << "  ys=y(it_,:);" << endl;
+  mDynamicModelFile << "  ys=y;" << endl
+                    << "  ll_index=M_.lead_lag_incidence(M_.maximum_endo_lag+1, :);" << endl;
 
   for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
     {
@@ -1323,7 +1323,9 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
         case BlockSimulationType::evaluateForward:
         case BlockSimulationType::evaluateBackward:
           mDynamicModelFile << "  [y, T, dr(" << blk + 1 << ").g1, dr(" << blk + 1 << ").g1_x, dr(" << blk + 1 << ").g1_xd, dr(" << blk + 1 << ").g1_o]=" << basename << ".block.dynamic_" << blk + 1 << "(y, x, params, steady_state, T, it_, true);" << endl
-                            << "  residual(y_index_eq)=ys(y_index)-y(it_, y_index);" << endl;
+                            << "  r=NaN(M_.endo_nbr,1);" << endl
+                            << "  r(nonzeros(ll_index))=ys(find(ll_index))-y(find(ll_index));" << endl
+                            << "  residual(y_index_eq)=r(y_index);" << endl;
           break;
         case BlockSimulationType::solveForwardSimple:
         case BlockSimulationType::solveBackwardSimple:
@@ -5496,7 +5498,7 @@ DynamicModel::isChecksumMatching(const string &basename, bool block) const
   // Write equation tags
   equation_tags.writeCheckSumInfo(buffer);
 
-  ExprNodeOutputType buffer_type = block ? ExprNodeOutputType::matlabDynamicBlockModel : ExprNodeOutputType::CDynamicModel;
+  ExprNodeOutputType buffer_type = ExprNodeOutputType::CDynamicModel;
 
   deriv_node_temp_terms_t tef_terms;
   temporary_terms_t temp_term_union;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index ac0f8dac..c9310a37 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -100,10 +100,8 @@ ExprNode::checkIfTemporaryTermThenWrite(ostream &output, ExprNodeOutputType outp
   // It is the responsibility of the caller to ensure that all temporary terms have their index
   assert(it2 != temporary_terms_idxs.end());
   output << "T" << LEFT_ARRAY_SUBSCRIPT(output_type)
-         << it2->second + ARRAY_SUBSCRIPT_OFFSET(output_type);
-  if (output_type == ExprNodeOutputType::matlabDynamicBlockModel)
-    output << ",it_";
-  output << RIGHT_ARRAY_SUBSCRIPT(output_type);
+         << it2->second + ARRAY_SUBSCRIPT_OFFSET(output_type)
+         << RIGHT_ARRAY_SUBSCRIPT(output_type);
 
   return true;
 }
@@ -1006,15 +1004,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
           output <<  "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
-        case ExprNodeOutputType::matlabDynamicBlockModel:
-          i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
-          if (lag > 0)
-            output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
-          else if (lag < 0)
-            output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_" << lag << ", " << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
-          else
-            output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_, " << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
-          break;
         case ExprNodeOutputType::matlabOutsideModel:
           output << "oo_.steady_state(" << tsid + 1 << ")";
           break;
@@ -1053,7 +1042,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         {
         case ExprNodeOutputType::juliaDynamicModel:
         case ExprNodeOutputType::matlabDynamicModel:
-        case ExprNodeOutputType::matlabDynamicBlockModel:
           if (lag > 0)
             output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i
                    << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -1112,7 +1100,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         {
         case ExprNodeOutputType::juliaDynamicModel:
         case ExprNodeOutputType::matlabDynamicModel:
-        case ExprNodeOutputType::matlabDynamicBlockModel:
           if (lag > 0)
             output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i
                    << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -2735,7 +2722,6 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       switch (output_type)
         {
         case ExprNodeOutputType::matlabDynamicModel:
-        case ExprNodeOutputType::matlabDynamicBlockModel:
           new_output_type = ExprNodeOutputType::matlabDynamicSteadyStateOperator;
           break;
         case ExprNodeOutputType::latexDynamicModel:
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 54d8b5c6..e6c8f1e3 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -83,7 +83,6 @@ enum class ExprNodeOutputType
   {
    matlabStaticModel, //!< Matlab code, static model
    matlabDynamicModel, //!< Matlab code, dynamic model
-   matlabDynamicBlockModel, //!< Matlab code, dynamic block-decomposed model
    CDynamicModel, //!< C code, dynamic model
    CStaticModel, //!< C code, static model
    juliaStaticModel, //!< Julia code, static model
@@ -107,7 +106,6 @@ isMatlabOutput(ExprNodeOutputType output_type)
   return output_type == ExprNodeOutputType::matlabStaticModel
     || output_type == ExprNodeOutputType::matlabDynamicModel
     || output_type == ExprNodeOutputType::matlabOutsideModel
-    || output_type == ExprNodeOutputType::matlabDynamicBlockModel
     || output_type == ExprNodeOutputType::matlabDynamicSteadyStateOperator
     || output_type == ExprNodeOutputType::steadyStateFile
     || output_type == ExprNodeOutputType::matlabDseries
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index bbe6f14d..8cb988cf 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -339,7 +339,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context) const
           catch (ExprNode::EvalException &e)
             {
               cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " (line " << equations_lineno[eq] << ") and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
-              d1->writeOutput(cerr, ExprNodeOutputType::matlabDynamicBlockModel, {}, {});
+              d1->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModel, {}, {});
               cerr << endl;
               exit(EXIT_FAILURE);
             }