diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 9ed5a3ba468483100b1b93982e95c047f9e6a78a..fe6bdda85a941a25352afd1450a87469aec9e274 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -488,53 +488,6 @@ DynamicModel::writeDynamicMCompatFile(const string &basename) const
   output.close();
 }
 
-void
-DynamicModel::writeDynamicJacobianNonZeroEltsFile(const string &basename) const
-{
-  vector<pair<int, int>> nzij_pred, nzij_current, nzij_fwrd; // pairs (tsid, equation)
-  for (const auto &[indices, d1] : derivatives[1])
-    {
-      if (getTypeByDerivID(indices[1]) != SymbolType::endogenous)
-        continue;
-      int tsid { getTypeSpecificIDByDerivID(indices[1]) };
-      int lag = getLagByDerivID(indices[1]);
-      if (lag == -1)
-        nzij_pred.emplace_back(tsid, indices[0]);
-      else if (lag == 0)
-        nzij_current.emplace_back(tsid, indices[0]);
-      else
-        nzij_fwrd.emplace_back(tsid, indices[0]);
-    }
-  sort(nzij_pred.begin(), nzij_pred.end());
-  sort(nzij_current.begin(), nzij_current.end());
-  sort(nzij_fwrd.begin(), nzij_fwrd.end());
-
-  const filesystem::path filename {packageDir(basename) / "dynamic_g1_nz.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 [nzij_pred, nzij_current, nzij_fwrd] = dynamic_g1_nz()" << endl
-         << "% Returns the coordinates of non-zero elements in the Jacobian, in column-major order, for each lead/lag (only for endogenous)" << endl;
-  auto print_nzij = [&output](const vector<pair<int, int>> &nzij, const string &name) {
-                      output << "  " << name << " = zeros(" << nzij.size() << ", 2, 'int32');" << endl;
-                      for (int idx{1};
-                           const auto &it : nzij)
-                        {
-                          output << "  " << name << "(" << idx << ",1)=" << it.second+1 << ';'
-                                 << " " << name << "(" << idx << ",2)=" << it.first+1 << ';' << endl;
-                          idx++;
-                        }
-                    };
-  print_nzij(nzij_pred, "nzij_pred");
-  print_nzij(nzij_current, "nzij_current");
-  print_nzij(nzij_fwrd, "nzij_fwrd");
-  output << "end" << endl;
-  output.close();
-}
-
 vector<pair<string, string>>
 DynamicModel::parseIncludeExcludeEquations(const string &inc_exc_option_value, bool exclude_eqs)
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 6d3dedfffeb472d9693b7fee01a4748f23ed41a3..a1f4b98d3cfeccb339b70960af81a7d687b5753f 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -373,10 +373,6 @@ public:
   template<bool julia>
   void writeParamsDerivativesFile(const string &basename) const;
 
-  //! Writes file containing coordinates of non-zero elements in the Jacobian
-  /*! Used by the perfect_foresight_problem MEX */
-  void writeDynamicJacobianNonZeroEltsFile(const string &basename) const;
-
   //! Creates mapping for variables and equations they are present in
   void createVariableMapping();
 
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 66bd8ba325a45dfcaa49e8a9ab813835a52105df..9c34b4702c0f7a48f5107a0536438872e26a19df 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -1050,8 +1050,6 @@ ModFile::writeMOutput(const string &basename, bool clear_all, bool clear_global,
           dynamic_model.writeDynamicFile(basename, use_dll, mexext, matlabroot, dynareroot, false);
 
           dynamic_model.writeParamsDerivativesFile<false>(basename);
-
-          dynamic_model.writeDynamicJacobianNonZeroEltsFile(basename);
         }
 
       // Create steady state file