diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index b9adcc219422f50479c519cc42373d81109f5f13..08f4f81882223aa9d43217c0b1ad97d55aaeddc7 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::matlabDynamicModelSparse;
+  constexpr ExprNodeOutputType local_output_type = ExprNodeOutputType::matlabDynamicBlockModel;
 
   for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
     {
@@ -254,15 +254,15 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
              << "%" << endl;
       if (simulation_type == BlockSimulationType::evaluateBackward
           || simulation_type == BlockSimulationType::evaluateForward)
-        output << "function [y, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)" << endl;
+        output << "function [y, T, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, jacobian_eval, y_kmin, periods)" << endl;
       else if (simulation_type == BlockSimulationType::solveForwardComplete
                || simulation_type == BlockSimulationType::solveBackwardComplete)
-        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
+        output << "function [residual, y, T, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, jacobian_eval)" << endl;
       else if (simulation_type == BlockSimulationType::solveBackwardSimple
                || simulation_type == BlockSimulationType::solveForwardSimple)
-        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
+        output << "function [residual, y, T, g1, g2, g3, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, jacobian_eval)" << endl;
       else
-        output << "function [residual, y, g1, g2, g3, b, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)" << endl;
+        output << "function [residual, y, T, g1, g2, g3, b, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, periods, jacobian_eval, y_kmin, y_size, Periods)" << endl;
 
       output << "  % ////////////////////////////////////////////////////////////////////////" << endl
              << "  % //" << string("                     Block ").substr(static_cast<int>(log10(blk + 1))) << blk+1
@@ -303,38 +303,6 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
       output << "  g2=0;" << endl
              << "  g3=0;" << endl;
 
-      // Declare global temp terms from this block and the previous ones
-      bool global_keyword_written = false;
-      for (int blk2 = 0; blk2 <= blk; blk2++)
-        for (auto &eq_tt : blocks_temporary_terms[blk2])
-          for (auto tt : eq_tt)
-            {
-              if (!global_keyword_written)
-                {
-                  output << "  global";
-                  global_keyword_written = true;
-                }
-              output << " ";
-              // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-              tt->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, eq_tt, {});
-            }
-      if (global_keyword_written)
-        output << ";" << endl;
-
-      // Initialize temp terms of this block
-      if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
-          || simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
-        {
-          output << "  " << "T_zeros = zeros(y_kmin+periods, 1);" << endl;
-          for (auto &eq_tt : blocks_temporary_terms[blk])
-            for (auto tt : eq_tt)
-              {
-                output << "  ";
-                // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-                tt->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, eq_tt, {});
-                output << " = T_zeros;" << endl;
-              }
-        }
       if (simulation_type == BlockSimulationType::solveBackwardSimple
           || simulation_type == BlockSimulationType::solveForwardSimple
           || simulation_type == BlockSimulationType::solveBackwardComplete
@@ -370,12 +338,12 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
                            for (auto it : blocks_temporary_terms[blk][eq])
                              {
                                if (dynamic_cast<AbstractExternalFunctionNode *>(it))
-                                 it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+                                 it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
 
                                output << indent_prefix;
-                               it->writeOutput(output, local_output_type, blocks_temporary_terms[blk][eq], {}, tef_terms);
+                               it->writeOutput(output, local_output_type, blocks_temporary_terms[blk][eq], blocks_temporary_terms_idxs, tef_terms);
                                output << " = ";
-                               it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms);
+                               it->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs, tef_terms);
                                temporary_terms.insert(it);
                                output << ";" << endl;
                              }
@@ -406,9 +374,9 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
                   exit(EXIT_FAILURE);
                 }
               output << indent_prefix;
-              lhs->writeOutput(output, local_output_type, temporary_terms, {});
+              lhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
               output << " = ";
-              rhs->writeOutput(output, local_output_type, temporary_terms, {});
+              rhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
               output << ";" << endl;
               break;
             case BlockSimulationType::solveBackwardSimple:
@@ -428,9 +396,9 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
               goto end;
             default:
             end:
-              lhs->writeOutput(output, local_output_type, temporary_terms, {});
+              lhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
               output << ") - (";
-              rhs->writeOutput(output, local_output_type, temporary_terms, {});
+              rhs->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
               output << ");" << endl;
             }
         }
@@ -448,7 +416,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
           auto [eq, var, lag] = indices;
           int jacob_col = blocks_jacob_cols_endo[blk].at({ var, lag });
           output << indent_prefix << "g1(" << eq+1 << ", " << jacob_col+1 << ") = ";
-          d->writeOutput(output, local_output_type, temporary_terms, {});
+          d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
           output << ";" << endl;
         }
       for (const auto &[indices, d] : blocks_derivatives_exo[blk])
@@ -456,7 +424,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
           auto [eq, var, lag] = indices;
           int jacob_col = blocks_jacob_cols_exo[blk].at({ var, lag });
           output << indent_prefix << "g1_x(" << eq+1 << ", " << jacob_col+1 << ") = ";
-          d->writeOutput(output, local_output_type, temporary_terms, {});
+          d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
           output << ";" << endl;
         }
       for (const auto &[indices, d] : blocks_derivatives_exo_det[blk])
@@ -464,7 +432,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
           auto [eq, var, lag] = indices;
           int jacob_col = blocks_jacob_cols_exo_det[blk].at({ var, lag });
           output << indent_prefix << "g1_xd(" << eq+1 << ", " << jacob_col+1 << ") = ";
-          d->writeOutput(output, local_output_type, temporary_terms, {});
+          d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
           output << ";" << endl;
         }
       for (const auto &[indices, d] : blocks_derivatives_other_endo[blk])
@@ -472,7 +440,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
           auto [eq, var, lag] = indices;
           int jacob_col = blocks_jacob_cols_other_endo[blk].at({ var, lag });
           output << indent_prefix << "g1_o(" << eq+1 << ", " << jacob_col+1 << ") = ";
-          d->writeOutput(output, local_output_type, temporary_terms, {});
+          d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
           output << ";" << endl;
         }
       output << indent_prefix << "varargout{1}=g1_x;" << endl
@@ -497,7 +465,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
               if (lag == 0)
                 {
                   output << "    g1(" << eq+1 << ", " << var+1-block_recursive_size << ") = ";
-                  d->writeOutput(output, local_output_type, temporary_terms, {});
+                  d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
                   output << ";" << endl;
                 }
             }
@@ -542,7 +510,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
                   else if (lag < 0)
                     output << "     g1(" << eq+1-block_recursive_size << "+Per_J_, "
                            << var+1-block_recursive_size << "+y_size*(it_" << lag-1 << ")) = ";
-                  d->writeOutput(output, local_output_type, temporary_terms, {});
+                  d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
                   output << ";" << endl;
                 }
             }
@@ -1446,30 +1414,10 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
 
   mDynamicModelFile << "function [varargout] = dynamic(options_, M_, oo_, varargin)" << endl
                     << "  g2=[];g3=[];" << endl;
-  //Temporary variables declaration
-  if (blocks_temporary_terms_idxs.size() > 0)
-    {
-      temporary_terms_t tt_all;
-      for (auto [tt, idx] : blocks_temporary_terms_idxs)
-        tt_all.insert(tt);
 
-      mDynamicModelFile << "  global";
-      for (auto tt : tt_all)
-        {
-          mDynamicModelFile << " ";
-          // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-          tt->writeOutput(mDynamicModelFile, ExprNodeOutputType::matlabStaticModelSparse, tt_all, {});
-        }
-      mDynamicModelFile << ";" << endl
-                        << "  T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);" << endl;
-      for (auto tt : tt_all)
-        {
-          mDynamicModelFile << "  ";
-          // In the following, "Static" is used to avoid getting the "(it_)" subscripting
-          tt->writeOutput(mDynamicModelFile, ExprNodeOutputType::matlabStaticModelSparse, tt_all, {});
-          mDynamicModelFile << "=T_init;" << 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 << "  y_kmin=M_.maximum_lag;" << endl
                     << "  y_kmax=M_.maximum_lead;" << endl
@@ -1519,22 +1467,22 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
         {
         case BlockSimulationType::evaluateForward:
         case BlockSimulationType::evaluateBackward:
-          mDynamicModelFile << "    [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 1, it_-1, 1);" << endl
+          mDynamicModelFile << "    [y, T, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, T, true, it_-1, 1);" << endl
                             << "    residual(y_index_eq)=ys(y_index)-y(it_, y_index);" << endl;
           break;
         case BlockSimulationType::solveForwardSimple:
         case BlockSimulationType::solveBackwardSimple:
-          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_, 1);" << endl
+          mDynamicModelFile << "    [r, y, T, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, T, it_, true);" << endl
                             << "    residual(y_index_eq)=r;" << endl;
           break;
         case BlockSimulationType::solveForwardComplete:
         case BlockSimulationType::solveBackwardComplete:
-          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_, 1);" << endl
+          mDynamicModelFile << "    [r, y, T, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, T, it_, true);" << endl
                             << "    residual(y_index_eq)=r;" << endl;
           break;
         case BlockSimulationType::solveTwoBoundariesComplete:
         case BlockSimulationType::solveTwoBoundariesSimple:
-          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" <<  block + 1 << "(y, x, params, steady_state, it_-" << max_lag << ", 1, " << max_lag << ", " << block_recursive << "," << "options_.periods" << ");" << endl
+          mDynamicModelFile << "    [r, y, T, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" <<  block + 1 << "(y, x, params, steady_state, T, it_-" << max_lag << ", true, " << max_lag << ", " << block_recursive << "," << "options_.periods" << ");" << endl
                             << "    residual(y_index_eq)=r(:,M_.maximum_lag+1);" << endl;
           break;
         default:
@@ -1601,7 +1549,7 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
                             << "  oo_.deterministic_simulation.block(blck_num).error = 0;" << endl
                             << "  oo_.deterministic_simulation.block(blck_num).iterations = 0;" << endl
                             << "  g1=[];g2=[];g3=[];" << endl
-                            << "  y=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);" << endl
+                            << "  [y, T] = " << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, T, false, y_kmin, periods);" << endl
                             << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
                             << "  if any(isnan(tmp) | isinf(tmp))" << endl
                             << "    disp(['Inf or Nan value during the evaluation of block " << block <<"']);" << endl
@@ -1627,7 +1575,7 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
                             << "  oo_.deterministic_simulation.block(blck_num).error = 0;" << endl
                             << "  oo_.deterministic_simulation.block(blck_num).iterations = 0;" << endl
                             << "  g1=[];g2=[];g3=[];" << endl
-                            << "  " << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);" << endl
+                            << "  [~, T] = " << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, T, false, y_kmin, periods);" << endl
                             << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
                             << "  if any(isnan(tmp) | isinf(tmp))" << endl
                             << "    disp(['Inf or Nan value during the evaluation of block " << block <<"']);" << endl
@@ -1655,8 +1603,8 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
                             << "  else" << endl
                             << "    blck_num = 1;" << endl
                             << "  end;" << endl
-                            << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
-                            << ", y, x, params, steady_state, [], y_index, " << nze
+                            << "  [y, T] = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
+                            << ", y, x, params, steady_state, T, y_index, " << nze
                             << ", options_.periods, " << (blocks[block].linear ? "true" : "false")
                             << ", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, true, true, false, M_, options_, oo_);" << endl
                             << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
@@ -1687,8 +1635,8 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
                             << "  else" << endl
                             << "    blck_num = 1;" << endl
                             << "  end;" << endl
-                            << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
-                            <<", y, x, params, steady_state, [], y_index, " << nze
+                            << "  [y, T] = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
+                            <<", y, x, params, steady_state, T, y_index, " << nze
                             <<", options_.periods, " << (blocks[block].linear ? "true" : "false")
                             <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, true, true, false, M_, options_, oo_);" << endl
                             << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
@@ -1717,8 +1665,8 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const
                             << "  else" << endl
                             << "    blck_num = 1;" << endl
                             << "  end;" << endl
-                            << "  [y oo_] = solve_two_boundaries('" << basename << ".block.dynamic_" <<  block + 1 << "'"
-                            <<", y, x, params, steady_state, y_index, " << nze
+                            << "  [y, T, oo_] = solve_two_boundaries('" << basename << ".block.dynamic_" <<  block + 1 << "'"
+                            <<", y, x, params, steady_state, T, y_index, " << nze
                             <<", options_.periods, " << blocks[block].max_lag
                             <<", " << blocks[block].max_lead
                             <<", " << (blocks[block].linear ? "true" : "false")
@@ -5885,7 +5833,7 @@ DynamicModel::isChecksumMatching(const string &basename, bool block) const
   // Write equation tags
   equation_tags.writeCheckSumInfo(buffer);
 
-  ExprNodeOutputType buffer_type = block ? ExprNodeOutputType::matlabDynamicModelSparse : ExprNodeOutputType::CDynamicModel;
+  ExprNodeOutputType buffer_type = block ? ExprNodeOutputType::matlabDynamicBlockModel : 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 4248a127be6af039143efbd1d8d36ba06b179b70..c8090afbb9d11f5f1f527063c10558b3098c8447 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -96,20 +96,15 @@ ExprNode::checkIfTemporaryTermThenWrite(ostream &output, ExprNodeOutputType outp
   if (auto it = temporary_terms.find(const_cast<ExprNode *>(this)); it == temporary_terms.end())
     return false;
 
-  if (output_type == ExprNodeOutputType::matlabDynamicModelSparse)
-    output << "T" << idx << "(it_)";
-  else
-    if (output_type == ExprNodeOutputType::matlabStaticModelSparse)
-      output << "T" << idx;
-    else
-      {
-        auto it2 = temporary_terms_idxs.find(const_cast<ExprNode *>(this));
-        // 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)
-               << RIGHT_ARRAY_SUBSCRIPT(output_type);
-      }
+  auto it2 = temporary_terms_idxs.find(const_cast<ExprNode *>(this));
+  // 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);
+
   return true;
 }
 
@@ -1008,11 +1003,10 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case ExprNodeOutputType::CStaticModel:
         case ExprNodeOutputType::juliaStaticModel:
         case ExprNodeOutputType::matlabStaticModel:
-        case ExprNodeOutputType::matlabStaticModelSparse:
           i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
           output <<  "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
-        case ExprNodeOutputType::matlabDynamicModelSparse:
+        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);
@@ -1059,7 +1053,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         {
         case ExprNodeOutputType::juliaDynamicModel:
         case ExprNodeOutputType::matlabDynamicModel:
-        case ExprNodeOutputType::matlabDynamicModelSparse:
+        case ExprNodeOutputType::matlabDynamicBlockModel:
           if (lag > 0)
             output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i
                    << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -1081,7 +1075,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case ExprNodeOutputType::CStaticModel:
         case ExprNodeOutputType::juliaStaticModel:
         case ExprNodeOutputType::matlabStaticModel:
-        case ExprNodeOutputType::matlabStaticModelSparse:
           output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case ExprNodeOutputType::matlabOutsideModel:
@@ -1119,7 +1112,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         {
         case ExprNodeOutputType::juliaDynamicModel:
         case ExprNodeOutputType::matlabDynamicModel:
-        case ExprNodeOutputType::matlabDynamicModelSparse:
+        case ExprNodeOutputType::matlabDynamicBlockModel:
           if (lag > 0)
             output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i
                    << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -1141,7 +1134,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case ExprNodeOutputType::CStaticModel:
         case ExprNodeOutputType::juliaStaticModel:
         case ExprNodeOutputType::matlabStaticModel:
-        case ExprNodeOutputType::matlabStaticModelSparse:
           output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case ExprNodeOutputType::matlabOutsideModel:
@@ -2743,7 +2735,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       switch (output_type)
         {
         case ExprNodeOutputType::matlabDynamicModel:
-        case ExprNodeOutputType::matlabDynamicModelSparse:
+        case ExprNodeOutputType::matlabDynamicBlockModel:
           new_output_type = ExprNodeOutputType::matlabDynamicSteadyStateOperator;
           break;
         case ExprNodeOutputType::latexDynamicModel:
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 4e9d1067f914604f0232d4abc09334c08a8a43bd..54d8b5c60f05f61401c415b41089dea1606f06a9 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -83,8 +83,7 @@ enum class ExprNodeOutputType
   {
    matlabStaticModel, //!< Matlab code, static model
    matlabDynamicModel, //!< Matlab code, dynamic model
-   matlabStaticModelSparse, //!< Matlab code, static block decomposed model
-   matlabDynamicModelSparse, //!< Matlab code, dynamic block decomposed model
+   matlabDynamicBlockModel, //!< Matlab code, dynamic block-decomposed model
    CDynamicModel, //!< C code, dynamic model
    CStaticModel, //!< C code, static model
    juliaStaticModel, //!< Julia code, static model
@@ -108,8 +107,7 @@ isMatlabOutput(ExprNodeOutputType output_type)
   return output_type == ExprNodeOutputType::matlabStaticModel
     || output_type == ExprNodeOutputType::matlabDynamicModel
     || output_type == ExprNodeOutputType::matlabOutsideModel
-    || output_type == ExprNodeOutputType::matlabStaticModelSparse
-    || output_type == ExprNodeOutputType::matlabDynamicModelSparse
+    || 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 f7616518fc4174411e204811939dd31efc3e5bef..c64e0330079c08edf5a6d3207687b741b6a1ea09 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::matlabDynamicModelSparse, {}, {});
+              d1->writeOutput(cerr, ExprNodeOutputType::matlabDynamicBlockModel, {}, {});
               cerr << endl;
               exit(EXIT_FAILURE);
             }
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 88404bd0163bc0c416d238cad26c633d3bb31e08..aa6be689412660844915f279eda58fd108b81f6e 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -152,7 +152,7 @@ StaticModel::writeStaticPerBlockMFiles(const string &basename) const
           || simulation_type == BlockSimulationType::evaluateForward)
         output << "function [y, T] = static_" << blk+1 << "(y, x, params, T)" << endl;
       else
-        output << "function [residual, y, g1, T] = static_" << blk+1 << "(y, x, params, T)" << endl;
+        output << "function [residual, y, T, g1] = static_" << blk+1 << "(y, x, params, T)" << endl;
 
       output << "  % ////////////////////////////////////////////////////////////////////////" << endl
              << "  % //" << string("                     Block ").substr(static_cast<int>(log10(blk + 1))) << blk+1
@@ -1730,7 +1730,7 @@ StaticModel::writeStaticBlockMFile(const string &basename) const
       exit(EXIT_FAILURE);
     }
 
-  output << "function [residual, g1, y, var_index, T] = static(nblock, y, x, params, T)" << endl
+  output << "function [residual, y, T, g1, var_index] = static(nblock, y, x, params, T)" << endl
          << "  residual = [];" << endl
          << "  g1 = [];" << endl
          << "  var_index = [];" << endl << endl
@@ -1756,7 +1756,7 @@ StaticModel::writeStaticBlockMFile(const string &basename) const
                  << "      y = y_tmp;" << endl;
         }
       else
-        output << "      [residual, y, g1, T] = " << basename << ".block.static_" << blk+1 << "(y, x, params, T);" << endl;
+        output << "      [residual, y, T, g1] = " << basename << ".block.static_" << blk+1 << "(y, x, params, T);" << endl;
 
     }
   output << "  end" << endl