diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 9ecb47a868389711289f8e485031086d12fc839e..1ae8342ea88e6f9860d8260665d8af0bfd765b96 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -20,6 +20,7 @@
 #include <sstream>
 #include <algorithm>
 #include <filesystem>
+#include <type_traits>
 
 #include "SparseMatrix.hh"
 
@@ -2229,8 +2230,8 @@ dynSparseMatrix::Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, d
 {
   size_t n = mxGetM(A_m);
   mxArray *z;
-  mxArray *rhs[2] = { A_m, b_m };
-  mexCallMATLAB(1, &z, 2, rhs, "mldivide");
+  mxArray *rhs[] = { A_m, b_m };
+  mexCallMATLAB(1, &z, std::extent_v<decltype(rhs)>, rhs, "mldivide");
   double *res = mxGetPr(z);
   if (is_two_boundaries)
     for (int i = 0; i < static_cast<int>(n); i++)
@@ -2468,20 +2469,20 @@ dynSparseMatrix::Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double
   size_t n = mxGetM(A_m);
   const char *field_names[] = {"droptol", "type"};
   mwSize dims[1] = { 1 };
-  mxArray *Setup = mxCreateStructArray(1, dims, 2, field_names);
+  mxArray *Setup = mxCreateStructArray(1, dims, std::extent_v<decltype(field_names)>, field_names);
   mxSetFieldByNumber(Setup, 0, 0, mxCreateDoubleScalar(lu_inc_tol));
   mxSetFieldByNumber(Setup, 0, 1, mxCreateString("ilutp"));
   mxArray *lhs0[2];
-  mxArray *rhs0[2] = { A_m, Setup };
-  if (mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"))
+  mxArray *rhs0[] = { A_m, Setup };
+  if (mexCallMATLAB(std::extent_v<decltype(lhs0)>, lhs0, std::extent_v<decltype(rhs0)>, rhs0, "ilu"))
     throw FatalException("In GMRES, the incomplete LU decomposition (ilu) has failed");
   mxArray *L1 = lhs0[0];
   mxArray *U1 = lhs0[1];
   /*[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);*/
-  mxArray *rhs[8] = { A_m, b_m, mxCreateDoubleScalar(Size), mxCreateDoubleScalar(1e-6),
-                      mxCreateDoubleScalar(static_cast<double>(n)), L1, U1, x0_m };
+  mxArray *rhs[] = { A_m, b_m, mxCreateDoubleScalar(Size), mxCreateDoubleScalar(1e-6),
+    mxCreateDoubleScalar(static_cast<double>(n)), L1, U1, x0_m };
   mxArray *lhs[2];
-  mexCallMATLAB(2, lhs, 8, rhs, "gmres");
+  mexCallMATLAB(std::extent_v<decltype(lhs)>, lhs, std::extent_v<decltype(rhs)>, rhs, "gmres");
   mxArray *z = lhs[0];
   mxArray *flag = lhs[1];
   double *flag1 = mxGetPr(flag);
@@ -2541,8 +2542,8 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
   if (preconditioner == 0)
     {
       mxArray *lhs0[1];
-      mxArray *rhs0[2] = { A_m, mxCreateDoubleScalar(0) };
-      mexCallMATLAB(1, lhs0, 2, rhs0, "spdiags");
+      mxArray *rhs0[] = { A_m, mxCreateDoubleScalar(0) };
+      mexCallMATLAB(std::extent_v<decltype(lhs0)>, lhs0, std::extent_v<decltype(rhs0)>, rhs0, "spdiags");
       mxArray *tmp = lhs0[0];
       double *tmp_val = mxGetPr(tmp);
       Diag = mxCreateSparse(n, n, n, mxREAL);
@@ -2563,15 +2564,15 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
       const char *field_names[] = {"type", "droptol", "milu", "udiag", "thresh"};
       const int type = 0, droptol = 1, milu = 2, udiag = 3, thresh = 4;
       mwSize dims[1] = { static_cast<mwSize>(1) };
-      mxArray *Setup = mxCreateStructArray(1, dims, 5, field_names);
+      mxArray *Setup = mxCreateStructArray(1, dims, std::extent_v<decltype(field_names)>, field_names);
       mxSetFieldByNumber(Setup, 0, type, mxCreateString("ilutp"));
       mxSetFieldByNumber(Setup, 0, droptol, mxCreateDoubleScalar(lu_inc_tol));
       mxSetFieldByNumber(Setup, 0, milu, mxCreateString("off"));
       mxSetFieldByNumber(Setup, 0, udiag, mxCreateDoubleScalar(0));
       mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(1));
       mxArray *lhs0[2];
-      mxArray *rhs0[2] = { A_m, Setup };
-      if (mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"))
+      mxArray *rhs0[] = { A_m, Setup };
+      if (mexCallMATLAB(std::extent_v<decltype(lhs0)>, lhs0, std::extent_v<decltype(rhs0)>, rhs0, "ilu"))
         throw FatalException{"In BiCGStab, the incomplete LU decomposition (ilu) has failed"};
       L1 = lhs0[0];
       U1 = lhs0[1];
@@ -2587,10 +2588,10 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
       for (int i = 0; i < static_cast<int>(n); i++)
         resid[i] = b[i] - resid[i];
       mxArray *lhs[1];
-      mxArray *rhs[2] = { L1, res };
-      mexCallMATLAB(1, lhs, 2, rhs, "mldivide");
-      mxArray *rhs2[2] = { U1, lhs[0] };
-      mexCallMATLAB(1, lhs, 2, rhs2, "mldivide");
+      mxArray *rhs[] = { L1, res };
+      mexCallMATLAB(std::extent_v<decltype(lhs)>, lhs, std::extent_v<decltype(rhs)>, rhs, "mldivide");
+      mxArray *rhs2[] = { U1, lhs[0] };
+      mexCallMATLAB(std::extent_v<decltype(lhs)>, lhs, std::extent_v<decltype(rhs2)>, rhs2, "mldivide");
       z = lhs[0];
       double *phat = mxGetPr(z);
       double *x0 = mxGetPr(x0_m);
@@ -2618,10 +2619,10 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
       if (preconditioner == 0)
         {
           /*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
-          mxArray *rhs[5] = { A_m, b_m, mxCreateDoubleScalar(1e-6),
-                              mxCreateDoubleScalar(static_cast<double>(n)), Diag };
+          mxArray *rhs[] = { A_m, b_m, mxCreateDoubleScalar(1e-6),
+            mxCreateDoubleScalar(static_cast<double>(n)), Diag };
           mxArray *lhs[2];
-          mexCallMATLAB(2, lhs, 5, rhs, "bicgstab");
+          mexCallMATLAB(std::extent_v<decltype(lhs)>, lhs, std::extent_v<decltype(rhs)>, rhs, "bicgstab");
           z = lhs[0];
           mxArray *flag = lhs[1];
           double *flag1 = mxGetPr(flag);
@@ -2634,10 +2635,10 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
       else if (preconditioner == 1)
         {
           /*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
-          mxArray *rhs[7] = { A_m, b_m, mxCreateDoubleScalar(1e-6),
-                              mxCreateDoubleScalar(static_cast<double>(n)), L1, U1, x0_m };
+          mxArray *rhs[] = { A_m, b_m, mxCreateDoubleScalar(1e-6),
+            mxCreateDoubleScalar(static_cast<double>(n)), L1, U1, x0_m };
           mxArray *lhs[2];
-          mexCallMATLAB(2, lhs, 7, rhs, "bicgstab");
+          mexCallMATLAB(std::extent_v<decltype(lhs)>, lhs, std::extent_v<decltype(rhs)>, rhs, "bicgstab");
           z = lhs[0];
           mxArray *flag = lhs[1];
           double *flag1 = mxGetPr(flag);
@@ -2693,7 +2694,7 @@ dynSparseMatrix::Singular_display(int block, int Size)
 {
   bool zero_solution;
   Simple_Init(Size, IM_i, zero_solution);
-  mxArray *rhs[1] = { mxCreateDoubleMatrix(Size, Size, mxREAL) };
+  mxArray *rhs[] = { mxCreateDoubleMatrix(Size, Size, mxREAL) };
   double *pind = mxGetPr(rhs[0]);
   for (int j = 0; j < Size * Size; j++)
     pind[j] = 0.0;
@@ -2710,7 +2711,7 @@ dynSparseMatrix::Singular_display(int block, int Size)
         }
     }
   mxArray *lhs[3];
-  mexCallMATLAB(3, lhs, 1, rhs, "svd");
+  mexCallMATLAB(std::extent_v<decltype(lhs)>, lhs, std::extent_v<decltype(rhs)>, rhs, "svd");
   mxArray *SVD_u = lhs[0];
   mxArray *SVD_s = lhs[1];
   double *SVD_ps = mxGetPr(SVD_s);
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index ab957f98dc6833b5da6688033b8aed87e1cd20a2..b17c553cef06635adf9560eb90ff6233c11c9558 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -20,6 +20,7 @@
 #include <ctime>
 #include <cmath>
 #include <cstring>
+#include <type_traits>
 
 #include "Interpreter.hh"
 #include "ErrorHandling.hh"
@@ -807,7 +808,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
                   jacob_exo_field_number = 1;
                   jacob_exo_det_field_number = 2;
                   mwSize dims[1] = { static_cast<mwSize>(nb_blocks) };
-                  plhs[1] = mxCreateStructArray(1, dims, 3, field_names);
+                  plhs[1] = mxCreateStructArray(1, dims, std::extent_v<decltype(field_names)>, field_names);
                 }
               else if (!mxIsStruct(block_structur))
                 {
diff --git a/mex/sources/k_order_perturbation/dynamic_m.cc b/mex/sources/k_order_perturbation/dynamic_m.cc
index 9c6035966fd7d9630fda8548205d30fb8098e02b..7846e441760e8c1d4ee8e9f0c854e751590466ee 100644
--- a/mex/sources/k_order_perturbation/dynamic_m.cc
+++ b/mex/sources/k_order_perturbation/dynamic_m.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2010-2022 Dynare Team
+ * Copyright © 2010-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -19,6 +19,7 @@
 
 #include <algorithm>
 #include <cassert>
+#include <type_traits>
 
 #include "dynare_exception.hh"
 
@@ -138,7 +139,7 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam
     std::string funcname = DynamicMFilename + "_g" + std::to_string(md.size()) + "_tt";
     mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m, steady_state_m, it_m };
 
-    int retVal = mexCallMATLAB(1, plhs, 6, prhs, funcname.c_str());
+    int retVal = mexCallMATLAB(std::extent_v<decltype(plhs)>, plhs, std::extent_v<decltype(prhs)>, prhs, funcname.c_str());
     if (retVal != 0)
       throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
 
@@ -151,7 +152,7 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam
     std::string funcname = DynamicMFilename + "_resid";
     mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m, steady_state_m, it_m, T_flag_m };
 
-    int retVal = mexCallMATLAB(1, plhs, 7, prhs, funcname.c_str());
+    int retVal = mexCallMATLAB(std::extent_v<decltype(plhs)>, plhs, std::extent_v<decltype(prhs)>, prhs, funcname.c_str());
     if (retVal != 0)
       throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
 
@@ -171,7 +172,7 @@ DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParam
       std::string funcname = DynamicMFilename + "_g" + std::to_string(i);
       mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m, steady_state_m, it_m, T_flag_m };
 
-      int retVal = mexCallMATLAB(1, plhs, 7, prhs, funcname.c_str());
+      int retVal = mexCallMATLAB(std::extent_v<decltype(plhs)>, plhs, std::extent_v<decltype(prhs)>, prhs, funcname.c_str());
       if (retVal != 0)
         throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
 
diff --git a/mex/sources/k_order_welfare/objective_m.cc b/mex/sources/k_order_welfare/objective_m.cc
index b3445753a375db9f59980a26dafee69dac276c43..2f0896a68698bff851b558b0f4c339c5f1a52392 100644
--- a/mex/sources/k_order_welfare/objective_m.cc
+++ b/mex/sources/k_order_welfare/objective_m.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2021 Dynare Team
+ * Copyright © 2021-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -19,6 +19,7 @@
 
 #include <algorithm>
 #include <cassert>
+#include <type_traits>
 
 #include "dynare_exception.hh"
 
@@ -92,7 +93,7 @@ ObjectiveMFile::eval(const Vector &y, const Vector &x, const Vector &modParams,
     std::string funcname = ObjectiveMFilename + "_g" + std::to_string(md.size()) + "_tt";
     mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m };
 
-    int retVal = mexCallMATLAB(1, plhs, 4, prhs, funcname.c_str());
+    int retVal = mexCallMATLAB(std::extent_v<decltype(plhs)>, plhs, std::extent_v<decltype(prhs)>, prhs, funcname.c_str());
     if (retVal != 0)
       throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
 
@@ -105,7 +106,7 @@ ObjectiveMFile::eval(const Vector &y, const Vector &x, const Vector &modParams,
     std::string funcname = ObjectiveMFilename + "_resid";
     mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m, T_flag_m };
 
-    int retVal = mexCallMATLAB(1, plhs, 5, prhs, funcname.c_str());
+    int retVal = mexCallMATLAB(std::extent_v<decltype(plhs)>, plhs, std::extent_v<decltype(prhs)>, prhs, funcname.c_str());
     if (retVal != 0)
       throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
 
@@ -119,7 +120,7 @@ ObjectiveMFile::eval(const Vector &y, const Vector &x, const Vector &modParams,
       std::string funcname = ObjectiveMFilename + "_g" + std::to_string(i);
       mxArray *plhs[1], *prhs[] = { T_m, y_m, x_m, params_m, T_flag_m };
 
-      int retVal = mexCallMATLAB(1, plhs, 5, prhs, funcname.c_str());
+      int retVal = mexCallMATLAB(std::extent_v<decltype(plhs)>, plhs, std::extent_v<decltype(prhs)>, prhs, funcname.c_str());
       if (retVal != 0)
         throw DynareException(__FILE__, __LINE__, "Trouble calling " + funcname);
 
diff --git a/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc b/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
index a3148239041b5a2fd3371aed03d5cd625bca7387..b5ad85e05f71860ec1c34246717592bda923a88f 100644
--- a/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
+++ b/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
@@ -189,7 +189,7 @@ DynamicModelMatlabCaller::eval(double *resid)
     std::string funcname {basename + ".sparse.dynamic_resid"};
     mxArray *plhs[3], *prhs[] = { y_mx, x_mx, params_mx, steady_state_mx };
 
-    mxArray *exception { mexCallMATLABWithTrap(3, plhs, 4, prhs, funcname.c_str()) };
+    mxArray *exception { mexCallMATLABWithTrap(std::extent_v<decltype(plhs)>, plhs, std::extent_v<decltype(prhs)>, prhs, funcname.c_str()) };
     if (exception)
       {
         error_msg = "An error occurred when calling " + funcname;
@@ -218,7 +218,7 @@ DynamicModelMatlabCaller::eval(double *resid)
       std::string funcname {basename + ".sparse.dynamic_g1"};
       mxArray *plhs[1], *prhs[] = { y_mx, x_mx, params_mx, steady_state_mx, g1_sparse_rowval_mx, g1_sparse_colval_mx, g1_sparse_colptr_mx, T_order_mx, T_mx };
 
-      mxArray *exception { mexCallMATLABWithTrap(1, plhs, 9, prhs, funcname.c_str()) };
+      mxArray *exception { mexCallMATLABWithTrap(std::extent_v<decltype(plhs)>, plhs, std::extent_v<decltype(prhs)>, prhs, funcname.c_str()) };
       if (exception)
         {
           error_msg = "An error occurred when calling " + funcname;