diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
index 3a79df7dbefa74b65c01c3acfd623926e4ae2974..fc18e31d469e599fefdfe985a956e04d2d1a5d24 100644
--- a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
+++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
@@ -20,6 +20,8 @@
 #ifndef _DYNAMICMODELAC_HH
 #define _DYNAMICMODELAC_HH
 
+#include <vector>
+
 #include "twod_matrix.hh"
 
 class DynamicModelAC
@@ -27,6 +29,6 @@ class DynamicModelAC
 public:
   virtual ~DynamicModelAC() = default;
   virtual void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
-                    Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) = 0;
+                    Vector &residual, std::vector<TwoDMatrix> &md) = 0;
 };
 #endif
diff --git a/mex/sources/k_order_perturbation/dynamic_dll.cc b/mex/sources/k_order_perturbation/dynamic_dll.cc
index ba8acd85c973d816e4e0ccfc6229066bc7949857..593e45a9ec348886c16993281578e7f30e63cdf1 100644
--- a/mex/sources/k_order_perturbation/dynamic_dll.cc
+++ b/mex/sources/k_order_perturbation/dynamic_dll.cc
@@ -94,21 +94,20 @@ DynamicModelDLL::~DynamicModelDLL()
 
 void
 DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady,
-                      Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false)
+                      Vector &residual, std::vector<TwoDMatrix> &md) noexcept(false)
 {
   dynamic_resid_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
   dynamic_resid(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), residual.base());
-  if (g1 || g2 || g3)
-    dynamic_g1_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
-  if (g1)
-    dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g1->base());
-  if (g2 || g3)
-    dynamic_g2_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
-  if (g2)
-    dynamic_g2(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g2->base());
-  if (g3)
+  dynamic_g1_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
+  dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), md[0].base());
+  if (md.size() >= 2)
+    {
+      dynamic_g2_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
+      dynamic_g2(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), md[1].base());
+    }
+  if (md.size() >= 3)
     {
       dynamic_g3_tt(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
-      dynamic_g3(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), g3->base());
+      dynamic_g3(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), md[2].base());
     }
 }
diff --git a/mex/sources/k_order_perturbation/dynamic_dll.hh b/mex/sources/k_order_perturbation/dynamic_dll.hh
index bfb1456dc69dd086b00de5ed3ec65196e3441efc..b29e5d66ea9aa037b8a12f8b170b0082ab92ad5e 100644
--- a/mex/sources/k_order_perturbation/dynamic_dll.hh
+++ b/mex/sources/k_order_perturbation/dynamic_dll.hh
@@ -66,6 +66,6 @@ public:
   virtual ~DynamicModelDLL();
 
   void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
-            Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) override;
+            Vector &residual, std::vector<TwoDMatrix> &md) override;
 };
 #endif
diff --git a/mex/sources/k_order_perturbation/dynamic_m.cc b/mex/sources/k_order_perturbation/dynamic_m.cc
index 0b05d53cfdb41b6d8dd63ecd31a5160798a61825..e129d169ebf5d98d76e9825bdaf1025f9c3cc09f 100644
--- a/mex/sources/k_order_perturbation/dynamic_m.cc
+++ b/mex/sources/k_order_perturbation/dynamic_m.cc
@@ -30,15 +30,15 @@ DynamicModelMFile::DynamicModelMFile(const std::string &modName) noexcept(false)
 }
 
 void
-DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm)
+DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix &tdm)
 {
   int totalCols = mxGetN(sparseMat);
   mwIndex *rowIdxVector = mxGetIr(sparseMat);
   mwSize sizeRowIdxVector = mxGetNzmax(sparseMat);
   mwIndex *colIdxVector = mxGetJc(sparseMat);
 
-  assert(tdm->ncols() == 3);
-  assert(tdm->nrows() == sizeRowIdxVector);
+  assert(tdm.ncols() == 3);
+  assert(tdm.nrows() == sizeRowIdxVector);
 
   double *ptr = mxGetPr(sparseMat);
 
@@ -48,9 +48,9 @@ DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat,
   for (int i = 0; i < totalCols; i++)
     for (int j = 0; j < static_cast<int>((colIdxVector[i+1]-colIdxVector[i])); j++, rind++)
       {
-        tdm->get(output_row, 0) = rowIdxVector[rind] + 1;
-        tdm->get(output_row, 1) = i + 1;
-        tdm->get(output_row, 2) = ptr[rind];
+        tdm.get(output_row, 0) = rowIdxVector[rind] + 1;
+        tdm.get(output_row, 1) = i + 1;
+        tdm.get(output_row, 2) = ptr[rind];
         output_row++;
       }
 
@@ -60,16 +60,16 @@ DynamicModelMFile::unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat,
      recognized as such by KordpDynare::populateDerivativesContainer() */
   while (output_row < static_cast<int>(sizeRowIdxVector))
     {
-      tdm->get(output_row, 0) = 0;
-      tdm->get(output_row, 1) = 0;
-      tdm->get(output_row, 2) = 0;
+      tdm.get(output_row, 0) = 0;
+      tdm.get(output_row, 1) = 0;
+      tdm.get(output_row, 2) = 0;
       output_row++;
     }
 }
 
 void
 DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady,
-                        Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false)
+                        Vector &residual, std::vector<TwoDMatrix> &md) noexcept(false)
 {
   constexpr int nlhs_dynamic = 4, nrhs_dynamic = 5;
   mxArray *prhs[nrhs_dynamic], *plhs[nlhs_dynamic];
@@ -91,14 +91,14 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam
 
   residual = Vector{plhs[0]};
 
-  assert(static_cast<int>(mxGetM(plhs[1])) == g1->nrows());
-  assert(static_cast<int>(mxGetN(plhs[1])) == g1->ncols());
-  std::copy_n(mxGetPr(plhs[1]), mxGetM(plhs[1])*mxGetN(plhs[1]), g1->base());
+  assert(static_cast<int>(mxGetM(plhs[1])) == md[0].nrows());
+  assert(static_cast<int>(mxGetN(plhs[1])) == md[0].ncols());
+  std::copy_n(mxGetPr(plhs[1]), mxGetM(plhs[1])*mxGetN(plhs[1]), md[0].base());
 
-  if (g2)
-    unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], g2);
-  if (g3)
-    unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[3], g3);
+  if (md.size() >= 2)
+    unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], md[1]);
+  if (md.size() >= 3)
+    unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[3], md[2]);
 
   for (int i = 0; i < nrhs_dynamic; i++)
     mxDestroyArray(prhs[i]);
diff --git a/mex/sources/k_order_perturbation/dynamic_m.hh b/mex/sources/k_order_perturbation/dynamic_m.hh
index 95d6dc3177b44ab2f8e3a5fd798b14f836bda019..0ae1e9aae3c6f9cc2feed523217e818141c7f932 100644
--- a/mex/sources/k_order_perturbation/dynamic_m.hh
+++ b/mex/sources/k_order_perturbation/dynamic_m.hh
@@ -33,11 +33,11 @@ class DynamicModelMFile : public DynamicModelAC
 {
 private:
   const std::string DynamicMFilename;
-  static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm);
+  static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix &tdm);
 public:
   explicit DynamicModelMFile(const std::string &modName);
   virtual ~DynamicModelMFile() = default;
   void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
-            Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) override;
+            Vector &residual, std::vector<TwoDMatrix> &md) override;
 };
 #endif
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc
index a4f5c5104addb7170f899cd15eb05a6a02a249e0..4523b0a85ef709e1f27f3ecffe80496b2eab4177 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.cc
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc
@@ -31,16 +31,13 @@ KordpDynare::KordpDynare(const std::vector<std::string> &endo,
                          int npred, int nforw, int nboth, const Vector &nnzd,
                          int nsteps, int norder,
                          Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
-                         const std::vector<int> &dr_order, const TwoDMatrix &llincidence,
-                         std::unique_ptr<TwoDMatrix> g1_arg, std::unique_ptr<TwoDMatrix> g2_arg,
-                         std::unique_ptr<TwoDMatrix> g3_arg) :
+                         const std::vector<int> &dr_order, const TwoDMatrix &llincidence) :
   nStat{nstat}, nBoth{nboth}, nPred{npred}, nForw{nforw}, nExog{nexog}, nPar{npar},
   nYs{npred + nboth}, nYss{nboth + nforw}, nY{nstat + npred + nboth + nforw},
   nJcols{nExog+nY+nYs+nYss}, NNZD{nnzd}, nSteps{nsteps},
   nOrder{norder}, journal{jr}, ySteady{ysteady}, params{inParams}, vCov{vcov},
   md{1}, dnl{*this, endo}, denl{*this, exo}, dsnl{*this, dnl, denl},
   ll_Incidence{llincidence},
-  g1p{std::move(g1_arg)}, g2p{std::move(g2_arg)}, g3p{std::move(g3_arg)},
   dynamicModelFile{std::move(dynamicModelFile_arg)}
 {
   computeJacobianPermutation(dr_order);
@@ -71,22 +68,16 @@ KordpDynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVect
 void
 KordpDynare::calcDerivativesAtSteady()
 {
-  if (!g1p)
+  if (dyn_md.empty())
     {
-      g1p = std::make_unique<TwoDMatrix>(nY, nJcols);
-      g1p->zeros();
+      dyn_md.emplace_back(nY, nJcols); // Allocate Jacobian
+      dyn_md.back().zeros();
 
-      if (nOrder > 1)
+      for (int i = 2; i <= nOrder; i++)
         {
-          // allocate space for sparse Hessian
-          g2p = std::make_unique<TwoDMatrix>(static_cast<int>(NNZD[1]), 3);
-          g2p->zeros();
-        }
-
-      if (nOrder > 2)
-        {
-          g3p = std::make_unique<TwoDMatrix>(static_cast<int>(NNZD[2]), 3);
-          g3p->zeros();
+          // Higher order derivatives, as sparse (3-column) matrices
+          dyn_md.emplace_back(static_cast<int>(NNZD[i-1]), 3);
+          dyn_md.back().zeros();
         }
 
       Vector xx(nexog());
@@ -97,21 +88,18 @@ KordpDynare::calcDerivativesAtSteady()
       Vector llxSteady(nJcols-nExog);
       LLxSteady(ySteady, llxSteady);
 
-      dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, g1p.get(), g2p.get(), g3p.get());
+      dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, dyn_md);
     }
 
-  populateDerivativesContainer(*g1p, 1);
-
-  if (nOrder > 1)
-    populateDerivativesContainer(*g2p, 2);
-
-  if (nOrder > 2)
-    populateDerivativesContainer(*g3p, 3);
+  for (int i = 1; i <= nOrder; i++)
+    populateDerivativesContainer(i);
 }
 
 void
-KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord)
+KordpDynare::populateDerivativesContainer(int ord)
 {
+  const TwoDMatrix &g = dyn_md[ord-1];
+
   // model derivatives FSSparseTensor instance
   auto mdTi = std::make_unique<FSSparseTensor>(ord, nJcols, nY);
 
@@ -249,6 +237,11 @@ KordpDynare::computeJacobianPermutation(const std::vector<int> &dr_order)
     dynToDynpp[dynppToDyn[i]] = i;
 }
 
+void
+KordpDynare::push_back_md(const mxArray *m)
+{
+  dyn_md.emplace_back(ConstTwoDMatrix{m});
+}
 
 DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector<std::string> names_arg)
   : names(std::move(names_arg))
@@ -263,3 +256,4 @@ DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const Dynare
   for (int i = 0; i < dynare.nexog(); i++)
     names.emplace_back(denl.getName(i));
 }
+
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.hh b/mex/sources/k_order_perturbation/k_ord_dynare.hh
index 87163d02ce9a04af2ecba183c61953f492e5fe4a..76722bcc7166a84242f135fe0ff4f038f0003307 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.hh
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.hh
@@ -95,15 +95,14 @@ private:
   Vector &ySteady;
   Vector &params;
   TwoDMatrix &vCov;
-  TensorContainer<FSSparseTensor> md; // ModelDerivatives
+  std::vector<TwoDMatrix> dyn_md; // Model derivatives, in Dynare form
+  TensorContainer<FSSparseTensor> md; // Model derivatives, in Dynare++ form
   DynareNameList dnl, denl;
   DynareStateNameList dsnl;
   const TwoDMatrix &ll_Incidence;
   std::vector<int> dynppToDyn; // Maps Dynare++ jacobian variable indices to Dynare ones
   std::vector<int> dynToDynpp; // Maps Dynare jacobian variable indices to Dynare++ ones
 
-  std::unique_ptr<TwoDMatrix> g1p, g2p, g3p;
-
   std::unique_ptr<DynamicModelAC> dynamicModelFile;
 public:
   KordpDynare(const std::vector<std::string> &endo,
@@ -112,9 +111,7 @@ public:
               int nforw, int nboth, const Vector &NNZD,
               int nSteps, int ord,
               Journal &jr, std::unique_ptr<DynamicModelAC> dynamicModelFile_arg,
-              const std::vector<int> &varOrder, const TwoDMatrix &ll_Incidence,
-              std::unique_ptr<TwoDMatrix> g1_arg,
-              std::unique_ptr<TwoDMatrix> g2_arg, std::unique_ptr<TwoDMatrix> g3_arg);
+              const std::vector<int> &varOrder, const TwoDMatrix &ll_Incidence);
 
   int
   nstat() const override
@@ -194,14 +191,16 @@ public:
     std::cerr << "KordpDynare::clone() not implemented" << std::endl;
     exit(EXIT_FAILURE);
   }
+  // Add model derivatives of a given order passed as argument to the MEX
+  void push_back_md(const mxArray *m);
 private:
   // Given the steady state in yS, returns in llxSteady the steady state extended with leads and lags
   void LLxSteady(const Vector &yS, Vector &llxSteady);
   /* Computes the permutations mapping back and forth between Dynare and
      Dynare++ orderings of variables */
   void computeJacobianPermutation(const std::vector<int> &var_order);
-  // Fills tensor container (with reordering) at a given order
-  void populateDerivativesContainer(const TwoDMatrix &g, int ord);
+  // Fills model derivatives in Dynare++ form (at a given order) given the Dynare form
+  void populateDerivativesContainer(int ord);
 };
 
 #endif
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 9f398db3c67772f146c39d278bbc27406b0652d8..b050b9b9298ab12f37bace41c04b89a5c40c7dd1 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -163,30 +163,6 @@ extern "C" {
     if (nEndo != static_cast<int>(endoNames.size()) || nExog != static_cast<int>(exoNames.size()))
       DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect size of M_.endo_names or M_.exo_names");
 
-    std::unique_ptr<TwoDMatrix> g1m, g2m, g3m;
-    if (nrhs > 3)
-      {
-        // Derivatives have been passed as arguments
-        const mxArray *g1 = prhs[3];
-        int m = static_cast<int>(mxGetM(g1));
-        int n = static_cast<int>(mxGetN(g1));
-        g1m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g1}});
-        if (nrhs > 4)
-          {
-            const mxArray *g2 = prhs[4];
-            int m = static_cast<int>(mxGetM(g2));
-            int n = static_cast<int>(mxGetN(g2));
-            g2m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g2}});
-            if (nrhs > 5)
-              {
-                const mxArray *g3 = prhs[5];
-                int m = static_cast<int>(mxGetM(g3));
-                int n = static_cast<int>(mxGetN(g3));
-                g3m = std::make_unique<TwoDMatrix>(m, n, Vector{ConstVector{g3}});
-              }
-          }
-      }
-
     const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
 
     try
@@ -206,11 +182,19 @@ extern "C" {
         KordpDynare dynare(endoNames, exoNames, nExog, nPar,
                            ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
                            NNZD, nSteps, kOrder, journal, std::move(dynamicModelFile),
-                           dr_order, llincidence,
-                           std::move(g1m), std::move(g2m), std::move(g3m));
+                           dr_order, llincidence);
 
-        // construct main K-order approximation class
+        // If model derivatives have been passed as arguments
+        if (nrhs > 3)
+          {
+            dynare.push_back_md(prhs[3]);
+            if (nrhs > 4)
+              dynare.push_back_md(prhs[4]);
+            if (nrhs > 5)
+              dynare.push_back_md(prhs[5]);
+          }
 
+        // construct main K-order approximation class
         Approximation app(dynare, journal, nSteps, false, qz_criterium);
         // run stochastic steady
         app.walkStochSteady();