diff --git a/dynare++/extern/matlab/dynare_simul.cc b/dynare++/extern/matlab/dynare_simul.cc
index 2271717cd78194590f61b3d145a2c99cf0c6d75b..2d84a511185821ffd113d514a74254df94145aaa 100644
--- a/dynare++/extern/matlab/dynare_simul.cc
+++ b/dynare++/extern/matlab/dynare_simul.cc
@@ -73,7 +73,7 @@ extern "C" {
     if (1 != ysteady_dim[1])
       DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of cols.\n");
 
-    mxArray *res = mxCreateDoubleMatrix(ny, nper, mxREAL);
+    plhs[1] = mxCreateDoubleMatrix(ny, nper, mxREAL);
 
     try
       {
@@ -84,16 +84,14 @@ extern "C" {
         UTensorPolynomial pol(ny, npred+nboth+nexog);
         for (int dim = 0; dim <= order; dim++)
           {
-            const mxArray *gk = prhs[11+dim];
-            const mwSize *const gk_dim = mxGetDimensions(gk);
+            ConstTwoDMatrix gk(prhs[11+dim]);
             FFSTensor ft(ny, npred+nboth+nexog, dim);
-            if (ft.ncols() != static_cast<int>(gk_dim[1]))
-              DYN_MEX_FUNC_ERR_MSG_TXT(("Wrong number of columns for folded tensor: got " + std::to_string(gk_dim[1]) + " but i want " + std::to_string(ft.ncols()) + '\n').c_str());
-            if (ft.nrows() != static_cast<int>(gk_dim[0]))
-              DYN_MEX_FUNC_ERR_MSG_TXT(("Wrong number of rows for folded tensor: got " + std::to_string(gk_dim[0]) + " but i want " + std::to_string(ft.nrows()) + '\n').c_str());
+            if (ft.ncols() != gk.ncols())
+              DYN_MEX_FUNC_ERR_MSG_TXT(("Wrong number of columns for folded tensor: got " + std::to_string(gk.ncols()) + " but i want " + std::to_string(ft.ncols()) + '\n').c_str());
+            if (ft.nrows() != gk.nrows())
+              DYN_MEX_FUNC_ERR_MSG_TXT(("Wrong number of rows for folded tensor: got " + std::to_string(gk.nrows()) + " but i want " + std::to_string(ft.nrows()) + '\n').c_str());
             ft.zeros();
-            ConstTwoDMatrix gk_mat(ft.nrows(), ft.ncols(), ConstVector{gk});
-            ft.add(1.0, gk_mat);
+            ft.add(1.0, gk);
             pol.insert(std::make_unique<UFSTensor>(ft));
           }
         // form the decision rule
@@ -106,9 +104,8 @@ extern "C" {
         // simulate and copy the results
         TwoDMatrix res_mat{dr.simulate(DecisionRule::emethod::horner, nper,
                                        ConstVector{ystart}, sr)};
-        TwoDMatrix res_tmp_mat{ny, nper, Vector{res}};
+        TwoDMatrix res_tmp_mat{plhs[1]};
         res_tmp_mat = const_cast<const TwoDMatrix &>(res_mat);
-        plhs[1] = res;
       }
     catch (const KordException &e)
       {
diff --git a/dynare++/sylv/cc/GeneralMatrix.hh b/dynare++/sylv/cc/GeneralMatrix.hh
index 715ccbdb7be742fbfe1c40ebc90534e5f49e8464..1c6e6857ddd90bb768dc39022f6e982ffce260ae 100644
--- a/dynare++/sylv/cc/GeneralMatrix.hh
+++ b/dynare++/sylv/cc/GeneralMatrix.hh
@@ -71,6 +71,12 @@ public:
   ConstGeneralMatrix(const GeneralMatrix &m, int i, int j, int nrows, int ncols);
   // Create submatrix (with data sharing)
   ConstGeneralMatrix(const ConstGeneralMatrix &m, int i, int j, int nrows, int ncols);
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+  explicit ConstGeneralMatrix(const mxArray *p)
+    : data(p), rows{static_cast<int>(mxGetM(p))}, cols{static_cast<int>(mxGetN(p))}, ld{rows}
+  {
+  }
+#endif
   virtual ~ConstGeneralMatrix() = default;
 
   ConstGeneralMatrix &operator=(const ConstGeneralMatrix &v) = delete;
@@ -208,6 +214,13 @@ public:
   // Create submatrix (with data sharing)
   GeneralMatrix(GeneralMatrix &m, int i, int j, int nrows, int ncols);
 
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+  explicit GeneralMatrix(mxArray *p)
+    : data(p), rows{static_cast<int>(mxGetM(p))}, cols{static_cast<int>(mxGetN(p))}, ld{rows}
+  {
+  }
+#endif
+
   virtual ~GeneralMatrix() = default;
   GeneralMatrix &operator=(const GeneralMatrix &m) = default;
   GeneralMatrix &operator=(GeneralMatrix &&m) = default;
diff --git a/dynare++/tl/cc/twod_matrix.hh b/dynare++/tl/cc/twod_matrix.hh
index 4241d51e0a620a2090e118743c4c331448fc826e..13ca1cb65fd99185ede1b1c9796c3e2c954a3b93 100644
--- a/dynare++/tl/cc/twod_matrix.hh
+++ b/dynare++/tl/cc/twod_matrix.hh
@@ -47,6 +47,11 @@ public:
     : ConstGeneralMatrix(m, first_row, first_col, rows, cols)
   {
   }
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+  explicit ConstTwoDMatrix(const mxArray *p) : ConstGeneralMatrix(p)
+  {
+  }
+#endif
   ~ConstTwoDMatrix() override = default;
 
   ConstTwoDMatrix &operator=(const ConstTwoDMatrix &v) = delete;
@@ -128,6 +133,13 @@ public:
     : GeneralMatrix(m, first_row, first_col, rows, cols)
   {
   }
+
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+  explicit TwoDMatrix(mxArray *p) : GeneralMatrix(p)
+  {
+  }
+#endif
+
   ~TwoDMatrix() override = default;
 
   TwoDMatrix &operator=(const TwoDMatrix &m) = default;
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 60f0829adcb78604e1d8672b56490bebdf9f1e88..9f398db3c67772f146c39d278bbc27406b0652d8 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -145,11 +145,9 @@ extern "C" {
                    [](double x) { return static_cast<int>(x)-1; });
 
     // the lag, current and lead blocks of the jacobian respectively
-    mxFldp = mxGetField(M_, 0, "lead_lag_incidence");
-    dim = static_cast<int>(mxGetN(mxFldp));
-    if (static_cast<int>(mxGetM(mxFldp)) != 3 || dim != nEndo)
+    TwoDMatrix llincidence(mxGetField(M_, 0, "lead_lag_incidence"));
+    if (llincidence.nrows() != 3 || llincidence.ncols() != nEndo)
       DYN_MEX_FUNC_ERR_MSG_TXT("Incorrect size of M_.lead_lag_incidence");
-    TwoDMatrix llincidence(3, dim, Vector{mxFldp});
 
     mxFldp = mxGetField(M_, 0, "NNZDerivatives");
     Vector NNZD{mxFldp};