diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 40b3d30cb13e1e1791dcddeb45a25b782cc75ecc..1ba99354df2e5b0437505ddf9745f32d050474a3 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3285,6 +3285,13 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool use_dll,
       create_directories(plusfolder);
       if (block && !use_dll)
         create_directories(plusfolder / "+block");
+
+      auto sparsefolder {plusfolder / "+sparse"};
+      create_directories(sparsefolder);
+      if (!use_dll)
+        create_directories(sparsefolder / "private");
+      if (block_decomposed)
+        create_directories(sparsefolder / "+block");
     }
   create_directories(model_dir / "bytecode");
 
@@ -3321,8 +3328,14 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool use_dll,
     }
 
   // Sparse representation
-  if (julia)
+  if (use_dll)
+    {
+      // TBD
+    }
+  else if (julia)
     writeSparseModelJuliaFiles<true>(basename);
+  else // MATLAB/Octave
+    writeSparseModelMFiles<true>(basename);
 
   writeSetAuxiliaryVariables(basename, julia);
 }
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 2936fdbd84a66cc5abf4e0fb45fed085b0c23ade..7d114f72bc24b607025d128a6c93522a13e55fde 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -118,9 +118,10 @@ private:
   //! Used for var_expectation and var_model
   map<string, set<int>> var_expectation_functions_to_write;
 
-  //! Writes dynamic model file (Matlab version)
+  // 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 version)
+  /* 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),
      then compiles it with the per-block functions into a single MEX */
@@ -131,7 +132,8 @@ private:
   // 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 version)
+  /* Writes the per-block dynamic files of block decomposed model (MATLAB/Octave
+     version, legacy representation) */
   void writeDynamicPerBlockMFiles(const string &basename) const;
   /* Writes and compiles the per-block dynamic files of block decomposed model
      (C version). Returns the list of paths to the compiled object files. */
@@ -206,9 +208,11 @@ private:
   //! Write reverse cross references
   void writeRevXrefs(ostream &output, const map<pair<int, int>, set<int>> &xrefmap, const string &type) const;
 
-  // Writes MATLAB/Octave wrapper function for computing residuals and derivatives at the same time
+  /* Writes MATLAB/Octave wrapper function for computing residuals and
+     derivatives at the same time (legacy representation) */
   void writeDynamicMWrapperFunction(const string &name, const string &ending) const;
-  // Helper for writing MATLAB/Octave functions for residuals/derivatives and their temporary terms
+  /* Helper for writing MATLAB/Octave functions for residuals/derivatives and
+     their temporary terms (legacy representation) */
   void writeDynamicMFileHelper(const string &basename,
                                const string &name, const string &retvalname,
                                const string &name_tt, size_t ttlen,
@@ -216,7 +220,8 @@ private:
                                const ostringstream &init_s, const ostringstream &end_s,
                                const ostringstream &s, const ostringstream &s_tt) const;
 
-  //! Create a legacy *_dynamic.m file for MATLAB/Octave not yet using the temporary terms array interface
+  /* Create the compatibility dynamic.m file for MATLAB/Octave not yet using
+     the temporary terms array interface (legacy representation) */
   void writeDynamicMCompatFile(const string &basename) const;
 
   //! Internal helper for the copy constructor and assignment operator
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 226b0e8d6e2e37f3758e67881db392bb1cc2d912..002a8f8978c23c5992a3b48cd558a3b5144f42df 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -375,6 +375,10 @@ protected:
   template<ExprNodeBytecodeOutputType output_type>
   void writeBytecodeModelEquations(BytecodeWriter &code_file, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms) const;
 
+  // Writes the sparse representation of the model in MATLAB/Octave
+  template<bool dynamic>
+  void writeSparseModelMFiles(const string &basename) const;
+
   // Writes the sparse representation of the model in Julia
   // Assumes that the directory <MODFILE>/model/julia/ already exists
   template<bool dynamic>
@@ -2445,3 +2449,164 @@ ModelTree::writeSparseModelJuliaFiles(const string &basename) const
     }
 }
 #endif
+
+template<bool dynamic>
+void
+ModelTree::writeSparseModelMFiles(const string &basename) const
+{
+  constexpr ExprNodeOutputType output_type {dynamic ? ExprNodeOutputType::matlabSparseDynamicModel : ExprNodeOutputType::matlabSparseStaticModel};
+  auto [d_sparse_output, tt_sparse_output] = writeModelFileHelper<output_type>();
+
+  const filesystem::path m_dir {packageDir(basename) / "+sparse"};
+  const filesystem::path private_m_dir {m_dir / "private"};
+  // TODO: when C++20 support is complete, mark the following strings constexpr
+  const string prefix { dynamic ? "dynamic_" : "static_" };
+  const string full_prefix { basename + ".sparse." + prefix };
+  const string ss_arg { dynamic ? ", steady_state" : "" };
+
+  size_t ttlen {temporary_terms_derivatives[0].size()};
+
+  ofstream output;
+  auto open_file = [&output](const filesystem::path &p)
+  {
+    output.open(p, ios::out | ios::binary);
+    if (!output.is_open())
+      {
+        cerr << "ERROR: Can't open file " << p.string() << " for writing" << endl;
+        exit(EXIT_FAILURE);
+      }
+  };
+
+  // Residuals (non-block)
+  open_file(private_m_dir / (prefix + "resid_tt.m"));
+  output << "function [T_order, T] = " << prefix << "resid_tt(y, x, params" << ss_arg << ", T_order, T)" << endl
+         << "if isempty(T_order)" << endl
+         << "    T_order = -1;" << endl
+         << "    T = NaN(" << ttlen << ", 1);" << endl
+         << "else if T_order >= 0" << endl
+         << "    return" << endl
+         << "end" << endl
+         << "T_order = 0;" << endl
+         << "if size(T, 1) < " << ttlen << endl
+         << "    T = resize(T, " << ttlen << ", 1);" << endl
+         << "end" << endl
+         << tt_sparse_output[0].str()
+         << "end" << endl;
+  output.close();
+
+  open_file(m_dir / (prefix + "resid.m"));
+  output << "function [residual, T_order, T] = " << prefix << "resid(y, x, params" << ss_arg << ", T_order, T)" << endl
+         << "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
+         << "residual = NaN(" << equations.size() << ", 1);" << endl
+         << d_sparse_output[0].str();
+  if constexpr(!dynamic)
+    output << "if ~isreal(residual)" << endl
+           << "    residual = real(residual)+imag(residual).^2;" << endl
+           << "end" << endl;
+  output << "end" << endl;
+  output.close();
+
+  // Jacobian (non-block)
+  ttlen += temporary_terms_derivatives[1].size();
+
+  open_file(private_m_dir / (prefix + "g1_tt.m"));
+  output << "function [T_order, T] = " << prefix << "g1_tt(y, x, params" << ss_arg << ", T_order, T)" << endl
+         << "if isempty(T_order)" << endl
+         << "    T_order = -1;" << endl
+         << "    T = NaN(" << ttlen << ", 1);" << endl
+         << "else if T_order >= 1" << endl
+         << "    return" << endl
+         << "end" << endl
+         << "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
+         << "T_order = 1;" << endl
+         << "if size(T, 1) < " << ttlen << endl
+         << "    T = resize(T, " << ttlen << ", 1);" << endl
+         << "end" << endl
+         << tt_sparse_output[1].str()
+         << "end" << endl;
+  output.close();
+
+  open_file(m_dir / (prefix + "g1.m"));
+  // NB: At first order, sparse indices are passed as extra arguments
+  output << "function [g1, T_order, T] = " << prefix << "g1(y, x, params" << ss_arg << ", sparse_rowval, sparse_colval, sparse_colptr, T_order, T)" << endl
+         << "[T_order, T] = " << full_prefix << "g1_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
+         << "g1_v = NaN(" << jacobian_sparse_column_major_order.size() << ", 1);" << endl
+         << d_sparse_output[1].str();
+  if constexpr(!dynamic)
+    output << "if ~isreal(g1_v)" << endl
+           << "    g1_v = real(g1_v)+2*imag(g1_v);" << endl
+           << "end" << endl;
+  output << "g1 = sparse(sparse_rowval, sparse_colval, g1_v, " << equations.size() << ", " << getJacobianColsNbr(true) << ");" << endl
+         << "end" << endl;
+  output.close();
+
+  // Higher-order derivatives (non-block)
+  for (int i {2}; i <= computed_derivs_order; i++)
+    {
+      ttlen += temporary_terms_derivatives[i].size();
+
+      open_file(private_m_dir / (prefix + "g" + to_string(i) + "_tt.m"));
+      output << "function T = " << prefix << "g" << i << "_tt(y, x, params" << ss_arg << ")" << endl
+             << "if isempty(T_order)" << endl
+             << "    T_order = -1;" << endl
+             << "    T = NaN(" << ttlen << ", 1);" << endl
+             << "else if T_order >= " << i << endl
+             << "    return" << endl
+             << "end" << endl
+             << "[T_order, T] = " << full_prefix << "g" << i-1 << "_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
+             << "T_order = " << i << ";" << endl
+             << "if size(T, 1) < " << ttlen << endl
+             << "    T = resize(T, " << ttlen << ", 1);" << endl
+             << "end" << endl
+             << tt_sparse_output[i].str()
+             << "end" << endl;
+      output.close();
+
+      open_file(m_dir / (prefix + "g" + to_string(i) + ".m"));
+      output << "function [g" << i << "_v, T_order, T] = " << prefix << "g" << i << "(y, x, params" << ss_arg << ", T_order, T)" << endl
+             << "[T_order, T] = " << full_prefix << "g" << i << "_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
+             << "g" << i << "_v = NaN(" << derivatives[i].size() << ", 1);" << endl
+             << d_sparse_output[i].str()
+             << "end" << endl;
+      output.close();
+    }
+
+  // Block decomposition
+  if (block_decomposed)
+    {
+      const filesystem::path block_dir {m_dir / "+block"};
+      temporary_terms_t temporary_terms_written;
+      for (int blk {0}; blk < static_cast<int>(blocks.size()); blk++)
+        {
+          const string funcname {prefix + to_string(blk+1)};
+          const BlockSimulationType simulation_type {blocks[blk].simulation_type};
+          const bool evaluate {simulation_type == BlockSimulationType::evaluateForward
+                               || simulation_type == BlockSimulationType::evaluateBackward};
+          const string resid_g1_arg {evaluate ? "" : ", residual, g1"};
+          open_file(block_dir / (funcname + ".m"));
+          output << "function [y, T" << resid_g1_arg << "] = " << funcname << "(y, x, params" << ss_arg << ", sparse_rowval, sparse_colval, sparse_colptr, T)" << endl;
+          if (!evaluate)
+            output << "residual=NaN(" << blocks[blk].mfs_size << ", 1);" << endl;
+
+          // Write residuals and temporary terms (incl. those for derivatives)
+          writePerBlockHelper<output_type>(blk, output, temporary_terms_written);
+
+          // Write Jacobian
+          if (!evaluate)
+            {
+              const bool one_boundary {simulation_type == BlockSimulationType::solveBackwardSimple
+                                       || simulation_type == BlockSimulationType::solveForwardSimple
+                                       || simulation_type == BlockSimulationType::solveBackwardComplete
+                                       || simulation_type == BlockSimulationType::solveForwardComplete};
+              output << "if nargout > 3" << endl
+                     << "    g1_v = NaN(" << blocks_jacobian_sparse_column_major_order[blk].size() << ", 1);" << endl;
+              writeSparsePerBlockJacobianHelper<output_type>(blk, output, temporary_terms_written);
+              output << "    g1 = sparse(sparse_rowval, sparse_colval, g1_v, " << blocks[blk].mfs_size << ", "
+                     << (one_boundary ? 1 : 3)*blocks[blk].mfs_size << ");" << endl
+                     << "end" << endl;
+            }
+          output << "end" << endl;
+          output.close();
+        }
+    }
+}
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index acb299d3ec7f4f41205b4291a96ef5d86c482b75..b1ace74507293ad067d8c61709dd2c1174d04cf0 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -628,6 +628,13 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool use_dll, c
       create_directories(plusfolder);
       if (block && !use_dll)
         create_directories(plusfolder / "+block");
+
+      auto sparsefolder {plusfolder / "+sparse"};
+      create_directories(sparsefolder);
+      if (!use_dll)
+        create_directories(sparsefolder / "private");
+      if (block_decomposed)
+        create_directories(sparsefolder / "+block");
     }
   create_directories(model_dir / "bytecode");
 
@@ -664,8 +671,14 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool use_dll, c
     }
 
   // Sparse representation
-  if (julia)
+  if (use_dll)
+    {
+      // TBD
+    }
+  else if (julia)
     writeSparseModelJuliaFiles<false>(basename);
+  else // MATLAB/Octave
+    writeSparseModelMFiles<false>(basename);
 
   writeSetAuxiliaryVariables(basename, julia);
 }
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index ae893c7aa9e1c7f0f893c97062043fa8eef663aa..36aadedb82d79ab6e3b03ad0b2f928988393acda 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -34,10 +34,11 @@ class DynamicModel;
 class StaticModel : public ModelTree
 {
 private:
-  //! Writes static model file (standard Matlab version)
+  // 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 version)
+  /* 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),
@@ -48,7 +49,8 @@ private:
   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 version)
+  /* Writes the per-block static files of block decomposed model (MATLAB/Octave
+     version, legacy representation) */
   void writeStaticPerBlockMFiles(const string &basename) const;
 
   /* Writes and compiles the per-block static files of block decomposed model
@@ -86,17 +88,20 @@ private:
 
   void computeChainRuleJacobian() override;
 
-  // Helper for writing MATLAB/Octave functions for residuals/derivatives and their temporary terms
+  /* Helper for writing MATLAB/Octave functions for residuals/derivatives and
+     their temporary terms (legacy representation) */
   void writeStaticMFileHelper(const string &basename,
                               const string &name, const string &retvalname,
                               const string &name_tt, size_t ttlen,
                               const string &previous_tt_name,
                               const ostringstream &init_s, const ostringstream &end_s,
                               const ostringstream &s, const ostringstream &s_tt) const;
-  // Writes MATLAB/Octave wrapper function for computing residuals and derivatives at the same time
+  /* Writes MATLAB/Octave wrapper function for computing residuals and
+     derivatives at the same time (legacy representation) */
   void writeStaticMWrapperFunction(const string &basename, const string &ending) const;
 
-  //! Create a legacy *_static.m file for MATLAB/Octave not yet using the temporary terms array interface
+  /* Create the compatibility static.m file for MATLAB/Octave not yet using the
+     temporary terms array interface (legacy representation) */
   void writeStaticMCompatFile(const string &name) const;
 
   int