diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
index e0da81e88..d4fd0f72f 100644
--- a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
+++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
@@ -33,40 +33,40 @@ DynamicModelAC::copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows,
       tdm->get(i, j) = dm[dmIdx++];
 }
 
-double *
-DynamicModelAC::unpackSparseMatrix(mxArray *sparseMat)
+void
+DynamicModelAC::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);
+
   double *ptr = mxGetPr(sparseMat);
-  double *newMat = (double *) malloc(sizeRowIdxVector*3*sizeof(double));
 
   int rind = 0;
-  int retvalind0 = 0;
-  int retvalind1 = sizeRowIdxVector;
-  int retvalind2 = sizeRowIdxVector*2;
+  int output_row = 0;
 
   for (int i = 0; i < totalCols; i++)
     for (int j = 0; j < (int) (colIdxVector[i+1]-colIdxVector[i]); j++, rind++)
       {
-        newMat[retvalind0++] = rowIdxVector[rind] + 1;
-        newMat[retvalind1++] = i + 1;
-        newMat[retvalind2++] = 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++;
       }
 
   /* If there are less elements than Nzmax (that might happen if some
      derivative is symbolically not zero but numerically zero at the evaluation
      point), then fill in the matrix with empty entries, that will be
      recognized as such by KordpDynare::populateDerivativesContainer() */
-  while (retvalind0 < (int) sizeRowIdxVector)
+  while (output_row < (int) sizeRowIdxVector)
     {
-      newMat[retvalind0++] = 0;
-      newMat[retvalind1++] = 0;
-      newMat[retvalind2++] = 0;
+      tdm->get(output_row, 0) = 0;
+      tdm->get(output_row, 1) = 0;
+      tdm->get(output_row, 2) = 0;
+      output_row++;
     }
-
-  return newMat;
 }
diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
index 0c5f63af6..e7f371d82 100644
--- a/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
+++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.hh
@@ -26,7 +26,7 @@ class DynamicModelAC
 {
 public:
   virtual ~DynamicModelAC() = default;
-  static double *unpackSparseMatrix(mxArray *sparseMatrix);
+  static void unpackSparseMatrixAndCopyIntoTwoDMatData(mxArray *sparseMat, TwoDMatrix *tdm);
   static void copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows, int cols);
   virtual void eval(const Vector &y, const Vector &x, const Vector &params, const Vector &ySteady,
                     Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) noexcept(false) = 0;
diff --git a/mex/sources/k_order_perturbation/dynamic_m.cc b/mex/sources/k_order_perturbation/dynamic_m.cc
index 0904dd4ae..162e98374 100644
--- a/mex/sources/k_order_perturbation/dynamic_m.cc
+++ b/mex/sources/k_order_perturbation/dynamic_m.cc
@@ -48,9 +48,9 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam
   residual = Vector(mxGetPr(plhs[0]), residual.skip(), (int) mxGetM(plhs[0]));
   copyDoubleIntoTwoDMatData(mxGetPr(plhs[1]), g1, (int) mxGetM(plhs[1]), (int) mxGetN(plhs[1]));
   if (g2 != nullptr)
-    copyDoubleIntoTwoDMatData(unpackSparseMatrix(plhs[2]), g2, (int) mxGetNzmax(plhs[2]), 3);
+    unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], g2);
   if (g3 != nullptr)
-    copyDoubleIntoTwoDMatData(unpackSparseMatrix(plhs[3]), g3, (int) mxGetNzmax(plhs[3]), 3);
+    unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[3], g3);
 
   for (int i = 0; i < nrhs_dynamic; i++)
     mxDestroyArray(prhs[i]);