diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
index dae58d34e9545ffab39833b246b0d58991b638c5..9630ed0f4f5edae9c3800699a96449c743fe106d 100644
--- a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
+++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc
@@ -37,40 +37,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 5449589252386a25abbd5b4ae56f5272f5569f00..c09d9d9807ee6ddb6c4648025dc92c3d7b075e2c 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();
-  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) throw (DynareException) = 0;
diff --git a/mex/sources/k_order_perturbation/dynamic_m.cc b/mex/sources/k_order_perturbation/dynamic_m.cc
index eb5fc69fd3723e7058c9749cb0ff19b830eaaaa4..3b617a1b0e3dd75a93c97a7d986b83243b48cf8c 100644
--- a/mex/sources/k_order_perturbation/dynamic_m.cc
+++ b/mex/sources/k_order_perturbation/dynamic_m.cc
@@ -52,9 +52,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 != NULL)
-    copyDoubleIntoTwoDMatData(unpackSparseMatrix(plhs[2]), g2, (int) mxGetNzmax(plhs[2]), 3);
+    unpackSparseMatrixAndCopyIntoTwoDMatData(plhs[2], g2);
   if (g3 != NULL)
-    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]);