diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index c53a4240edddc678508eda95f88f05355f942c61..f29d4283f7a54aed7d6dea248ecacf368d7f4cf9 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -2280,6 +2280,35 @@ DynamicModel::writeDynamicModelHelper(const string &name, const string &retvalna
   output.close();
 }
 
+void
+DynamicModel::writeDynamicMatlabCompatLayer(const string &name) const
+{
+  string filename = name + ".m";
+  ofstream output;
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  int ntt = temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size();
+
+  output << "function [residual, g1, g2, g3] = " << name << "(y, x, params, steady_state, it_)" << endl
+         << "    T = NaN(" << ntt << ", 1);" << endl
+         << "    if nargout <= 1" << endl
+         << "        residual = " << name << "_resid(T, y, x, params, steady_state, it_, true);" << endl
+         << "    elseif nargout == 2" << endl
+         << "        [residual, g1] = " << name << "_resid_g1(T, y, x, params, steady_state, it_, true);" << endl
+         << "    elseif nargout == 3" << endl
+         << "        [residual, g1, g2] = " << name << "_resid_g1_g2(T, y, x, params, steady_state, it_, true);" << endl
+         << "    else" << endl
+         << "        [residual, g1, g2, g3] = " << name << "_resid_g1_g2_g3(T, y, x, params, steady_state, it_, true);" << endl
+         << "    end" << endl
+         << "end" << endl;
+
+  output.close();
+}
+
 void
 DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const
 {
@@ -2572,6 +2601,8 @@ DynamicModel::writeDynamicModel(const string &dynamic_basename, ostream &Dynamic
                               init_output, end_output,
                               third_derivatives_output, third_derivatives_tt_output);
       writeWrapperFunctions(dynamic_basename, "g3");
+
+      writeDynamicMatlabCompatLayer(dynamic_basename);
     }
         /*
 
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index ce90586a1ae8f28ec7d548ce666adfff1120f210..28fd2c1a168586b01b31f3b1c0127ab8e529316c 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -245,6 +245,9 @@ private:
                                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
+  void writeDynamicMatlabCompatLayer(const string &name) const;
+
 public:
   DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_argx);
   //! Adds a variable node
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 7b50fd868310686a56b0b2420653a27ae6cce58e..e90b2ff4642db5da2e692f11395748337c32d16e 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1294,6 +1294,35 @@ StaticModel::writeStaticModelHelper(const string &name, const string &retvalname
   output.close();
 }
 
+void
+StaticModel::writeStaticMatlabCompatLayer(const string &name) const
+{
+  string filename = name + ".m";
+  ofstream output;
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+  int ntt = temporary_terms_mlv.size() + temporary_terms_res.size() + temporary_terms_g1.size() + temporary_terms_g2.size() + temporary_terms_g3.size();
+
+  output << "function [residual, g1, g2, g3] = " << name << "(y, x, params)" << endl
+         << "    T = NaN(" << ntt << ", 1);" << endl
+         << "    if nargout <= 1" << endl
+         << "        residual = " << name << "_resid(T, y, x, params, true);" << endl
+         << "    elseif nargout == 2" << endl
+         << "        [residual, g1] = " << name << "_resid_g1(T, y, x, params, true);" << endl
+         << "    elseif nargout == 3" << endl
+         << "        [residual, g1, g2] = " << name << "_resid_g1_g2(T, y, x, params, true);" << endl
+         << "    else" << endl
+         << "        [residual, g1, g2, g3] = " << name << "_resid_g1_g2_g3(T, y, x, params, true);" << endl
+         << "    end" << endl
+         << "end" << endl;
+
+  output.close();
+}
+
 void
 StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const
 {
@@ -1598,6 +1627,8 @@ StaticModel::writeStaticModel(const string &basename,
                              init_output, end_output,
                              third_derivatives_output, third_derivatives_tt_output);
       writeWrapperFunctions(static_name, "g3");
+
+      writeStaticMatlabCompatLayer(static_name);
     }
   else if (output_type == oCStaticModel)
     {
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 54394d0e353dad77392d95377d0d4082b6e43318..d229ce5380f0838649b276e6879eec8ffe0f303c 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -155,6 +155,10 @@ protected:
                               const ostringstream &init_s, const ostringstream &end_s,
                               const ostringstream &s, const ostringstream &s_tt) const;
   void writeWrapperFunctions(const string &basename, const string &ending) const;
+
+  //! Create a legacy *_static.m file for Matlab/Octave not yet using the temporary terms array interface
+  void writeStaticMatlabCompatLayer(const string &name) const;
+
   void writeStaticModel(ostream &DynamicOutput, bool use_dll, bool julia) const;
   void writeStaticModel(const string &dynamic_basename, bool use_dll, bool julia) const;
 public: