diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh
index dd462e003..cb274d5de 100644
--- a/mex/sources/bytecode/ErrorHandling.hh
+++ b/mex/sources/bytecode/ErrorHandling.hh
@@ -330,7 +330,7 @@ public:
       {
         nb_endo = 0;
         endo_name_length = 0;
-        P_endo_names = NULL;
+        P_endo_names = nullptr;
       }
     else
       {
@@ -342,7 +342,7 @@ public:
       {
         nb_exo = 0;
         exo_name_length = 0;
-        P_exo_names = NULL;
+        P_exo_names = nullptr;
       }
     else
       {
@@ -354,7 +354,7 @@ public:
       {
         nb_param = 0;
         param_name_length = 0;
-        P_param_names = NULL;
+        P_param_names = nullptr;
       }
     else
       {
@@ -602,7 +602,7 @@ public:
     double ll;
     ExpressionType equation_type = ExpressionType::TemporaryTerm;
     size_t found;
-    double *jacob = NULL, *jacob_other_endo = NULL, *jacob_exo = NULL, *jacob_exo_det = NULL;
+    double *jacob = nullptr, *jacob_other_endo = nullptr, *jacob_exo = nullptr, *jacob_exo_det = nullptr;
     ExternalFunctionType function_type = ExternalFunctionType::withoutDerivative;
 
     if (evaluate)
@@ -1856,9 +1856,9 @@ public:
                           }
                         input_arguments[nb_input_arguments+nb_add_input_arguments] = vv;
 #ifdef DEBUG
-                        mexCallMATLAB(0, NULL, 1, &input_arguments[0], "disp");
-                        mexCallMATLAB(0, NULL, 1, &input_arguments[1], "disp");
-                        mexCallMATLAB(0, NULL, 1, &input_arguments[2], "celldisp");
+                        mexCallMATLAB(0, nullptr, 1, &input_arguments[0], "disp");
+                        mexCallMATLAB(0, nullptr, 1, &input_arguments[1], "disp");
+                        mexCallMATLAB(0, nullptr, 1, &input_arguments[2], "celldisp");
                         mexPrintf("OK\n");
                         mexEvalString("drawnow;");
 #endif
@@ -1943,9 +1943,9 @@ public:
                           }
                         input_arguments[nb_input_arguments+nb_add_input_arguments] = vv;
 #ifdef DEBUG
-                        mexCallMATLAB(0, NULL, 1, &input_arguments[0], "disp");
-                        mexCallMATLAB(0, NULL, 1, &input_arguments[1], "disp");
-                        mexCallMATLAB(0, NULL, 1, &input_arguments[2], "celldisp");
+                        mexCallMATLAB(0, nullptr, 1, &input_arguments[0], "disp");
+                        mexCallMATLAB(0, nullptr, 1, &input_arguments[1], "disp");
+                        mexCallMATLAB(0, nullptr, 1, &input_arguments[2], "celldisp");
                         mexPrintf("OK\n");
                         mexEvalString("drawnow;");
 #endif
@@ -2224,7 +2224,7 @@ public:
   inline
   test_mxMalloc(void *z, int line, string file, string func, int amount)
   {
-    if (z == NULL && (amount > 0))
+    if (z == nullptr && (amount > 0))
       {
         ostringstream tmp;
         tmp << " mxMalloc: out of memory " << amount << " bytes required at line " << line << " in function " << func << " (file " << file;
diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc
index 9ed6313bf..4b9320397 100644
--- a/mex/sources/bytecode/Evaluate.cc
+++ b/mex/sources/bytecode/Evaluate.cc
@@ -117,7 +117,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
   bool go_on = true;
   double ll;
   double rr;
-  double *jacob = NULL, *jacob_other_endo = NULL, *jacob_exo = NULL, *jacob_exo_det = NULL;
+  double *jacob = nullptr, *jacob_other_endo = nullptr, *jacob_exo = nullptr, *jacob_exo_det = nullptr;
   EQN_block = block_num;
   stack<double> Stack;
   ExternalFunctionType function_type = ExternalFunctionType::withoutDerivative;
@@ -1232,9 +1232,9 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
                     }
                   input_arguments[nb_input_arguments+nb_add_input_arguments] = vv;
 #ifdef DEBUG
-                  mexCallMATLAB(0, NULL, 1, &input_arguments[0], "disp");
-                  mexCallMATLAB(0, NULL, 1, &input_arguments[1], "disp");
-                  mexCallMATLAB(0, NULL, 1, &input_arguments[2], "celldisp");
+                  mexCallMATLAB(0, nullptr, 1, &input_arguments[0], "disp");
+                  mexCallMATLAB(0, nullptr, 1, &input_arguments[1], "disp");
+                  mexCallMATLAB(0, nullptr, 1, &input_arguments[2], "celldisp");
                   mexPrintf("OK\n");
                   mexEvalString("drawnow;");
 #endif
@@ -1298,9 +1298,9 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int
                     }
                   input_arguments[nb_input_arguments+nb_add_input_arguments] = vv;
 #ifdef DEBUG
-                  mexCallMATLAB(0, NULL, 1, &input_arguments[0], "disp");
-                  mexCallMATLAB(0, NULL, 1, &input_arguments[1], "disp");
-                  mexCallMATLAB(0, NULL, 1, &input_arguments[2], "celldisp");
+                  mexCallMATLAB(0, nullptr, 1, &input_arguments[0], "disp");
+                  mexCallMATLAB(0, nullptr, 1, &input_arguments[1], "disp");
+                  mexCallMATLAB(0, nullptr, 1, &input_arguments[2], "celldisp");
                   mexPrintf("OK\n");
                   mexEvalString("drawnow;");
 #endif
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index c5dd79fb9..b34bcc066 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -65,7 +65,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
   y_decal = y_decal_arg;
   markowitz_c = markowitz_c_arg;
   filename = filename_arg;
-  T = NULL;
+  T = nullptr;
   minimal_solving_periods = minimal_solving_periods_arg;
   stack_solve_algo = stack_solve_algo_arg;
   solve_algo = solve_algo_arg;
@@ -775,7 +775,7 @@ Interpreter::MainLoop(string bin_basename, CodeLoad code, bool evaluate, int blo
             mxFree(T);
           if (global_temporary_terms)
             {
-              if (GlobalTemporaryTerms == NULL)
+              if (GlobalTemporaryTerms == nullptr)
                 {
                   mexPrintf("GlobalTemporaryTerms is NULL\n"); mexEvalString("drawnow;");
                 }
diff --git a/mex/sources/bytecode/Mem_Mngr.cc b/mex/sources/bytecode/Mem_Mngr.cc
index 9695b24bd..e706bd92c 100644
--- a/mex/sources/bytecode/Mem_Mngr.cc
+++ b/mex/sources/bytecode/Mem_Mngr.cc
@@ -41,8 +41,8 @@ Mem_Mngr::init_Mem()
   Chunk_Stack.clear();
   CHUNK_SIZE = 0;
   Nb_CHUNK = 0;
-  NZE_Mem = NULL;
-  NZE_Mem_add = NULL;
+  NZE_Mem = nullptr;
+  NZE_Mem_add = nullptr;
   CHUNK_heap_pos = 0;
   NZE_Mem_Allocated.clear();
 }
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index 0c013cb9b..54136beec 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -88,8 +88,8 @@ t_umfpack_dl_defaults umfpack_dl_defaults;
 
 dynSparseMatrix::dynSparseMatrix()
 {
-  pivotva = NULL;
-  g_save_op = NULL;
+  pivotva = nullptr;
+  g_save_op = nullptr;
   g_nop_all = 0;
   mem_mngr.init_Mem();
   symbolic = true;
@@ -103,8 +103,8 @@ dynSparseMatrix::dynSparseMatrix()
   restart = 0;
   IM_i.clear();
   lu_inc_tol = 1e-10;
-  Symbolic = NULL;
-  Numeric = NULL;
+  Symbolic = nullptr;
+  Numeric = nullptr;
 #ifdef _MSC_VER
   // Get a handle to the DLL module.
   hinstLib = LoadLibrary(TEXT("libmwumfpack.dll"));
@@ -187,8 +187,8 @@ dynSparseMatrix::dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, con
                                  ) :
   Evaluate(y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, periods_arg, minimal_solving_periods_arg, slowc_arg)
 {
-  pivotva = NULL;
-  g_save_op = NULL;
+  pivotva = nullptr;
+  g_save_op = nullptr;
   g_nop_all = 0;
   mem_mngr.init_Mem();
   symbolic = true;
@@ -202,8 +202,8 @@ dynSparseMatrix::dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, con
   restart = 0;
   IM_i.clear();
   lu_inc_tol = 1e-10;
-  Symbolic = NULL;
-  Numeric = NULL;
+  Symbolic = nullptr;
+  Numeric = nullptr;
 #ifdef CUDA
   CUDA_device = CUDA_device_arg;
   cublas_handle = cublas_handle_arg;
@@ -380,28 +380,28 @@ dynSparseMatrix::At_Col(int c, int lag, NonZeroElem **first)
 void
 dynSparseMatrix::Delete(const int r, const int c)
 {
-  NonZeroElem *first = FNZE_R[r], *firsta = NULL;
+  NonZeroElem *first = FNZE_R[r], *firsta = nullptr;
 
   while (first->c_index != c)
     {
       firsta = first;
       first = first->NZE_R_N;
     }
-  if (firsta != NULL)
+  if (firsta != nullptr)
     firsta->NZE_R_N = first->NZE_R_N;
   if (first == FNZE_R[r])
     FNZE_R[r] = first->NZE_R_N;
   NbNZRow[r]--;
 
   first = FNZE_C[c];
-  firsta = NULL;
+  firsta = nullptr;
   while (first->r_index != r)
     {
       firsta = first;
       first = first->NZE_C_N;
     }
 
-  if (firsta != NULL)
+  if (firsta != nullptr)
     firsta->NZE_C_N = first->NZE_C_N;
   if (first == FNZE_C[c])
     FNZE_C[c] = first->NZE_C_N;
@@ -465,7 +465,7 @@ dynSparseMatrix::Insert(const int r, const int c, const int u_index, const int l
   NonZeroElem *firstn, *first, *firsta, *a;
   firstn = mem_mngr.mxMalloc_NZE();
   first = FNZE_R[r];
-  firsta = NULL;
+  firsta = nullptr;
   while (first->c_index < c && (a = first->NZE_R_N))
     {
       firsta = first;
@@ -479,18 +479,18 @@ dynSparseMatrix::Insert(const int r, const int c, const int u_index, const int l
     {
       if (first == FNZE_R[r])
         FNZE_R[r] = firstn;
-      if (firsta != NULL)
+      if (firsta != nullptr)
         firsta->NZE_R_N = firstn;
       firstn->NZE_R_N = first;
     }
   else
     {
       first->NZE_R_N = firstn;
-      firstn->NZE_R_N = NULL;
+      firstn->NZE_R_N = nullptr;
     }
   NbNZRow[r]++;
   first = FNZE_C[c];
-  firsta = NULL;
+  firsta = nullptr;
   while (first->r_index < r && (a = first->NZE_C_N))
     {
       firsta = first;
@@ -500,14 +500,14 @@ dynSparseMatrix::Insert(const int r, const int c, const int u_index, const int l
     {
       if (first == FNZE_C[c])
         FNZE_C[c] = firstn;
-      if (firsta != NULL)
+      if (firsta != nullptr)
         firsta->NZE_C_N = firstn;
       firstn->NZE_C_N = first;
     }
   else
     {
       first->NZE_C_N = firstn;
-      firstn->NZE_C_N = NULL;
+      firstn->NZE_C_N = nullptr;
     }
 
   NbNZCol[c]++;
@@ -654,7 +654,7 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
   test_mxMalloc(line_done, __LINE__, __FILE__, __func__, Size*sizeof(bool));
 
   mem_mngr.init_CHUNK_BLCK_SIZE(u_count);
-  g_save_op = NULL;
+  g_save_op = nullptr;
   g_nop_all = 0;
   i = Size*sizeof(NonZeroElem *);
   FNZE_R = static_cast<NonZeroElem **>(mxMalloc(i));
@@ -675,10 +675,10 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
   for (i = 0; i < Size; i++)
     {
       line_done[i] = false;
-      FNZE_C[i] = NULL;
-      FNZE_R[i] = NULL;
-      temp_NZE_C[i] = 0;
-      temp_NZE_R[i] = 0;
+      FNZE_C[i] = nullptr;
+      FNZE_R[i] = nullptr;
+      temp_NZE_C[i] = nullptr;
+      temp_NZE_R[i] = nullptr;
       NbNZRow[i] = 0;
       NbNZCol[i] = 0;
     }
@@ -693,19 +693,19 @@ dynSparseMatrix::Simple_Init(int Size, map<pair<pair<int, int>, int>, int> &IM,
           NbNZRow[eq]++;
           NbNZCol[var]++;
           first = mem_mngr.mxMalloc_NZE();
-          first->NZE_C_N = NULL;
-          first->NZE_R_N = NULL;
+          first->NZE_C_N = nullptr;
+          first->NZE_R_N = nullptr;
           first->u_index = u_count1;
           first->r_index = eq;
           first->c_index = var;
           first->lag_index = lag;
-          if (FNZE_R[eq] == NULL)
+          if (FNZE_R[eq] == nullptr)
             FNZE_R[eq] = first;
-          if (FNZE_C[var] == NULL)
+          if (FNZE_C[var] == nullptr)
             FNZE_C[var] = first;
-          if (temp_NZE_R[eq] != NULL)
+          if (temp_NZE_R[eq] != nullptr)
             temp_NZE_R[eq]->NZE_R_N = first;
-          if (temp_NZE_C[var] != NULL)
+          if (temp_NZE_C[var] != nullptr)
             temp_NZE_C[var]->NZE_C_N = first;
           temp_NZE_R[eq] = first;
           temp_NZE_C[var] = first;
@@ -1065,7 +1065,7 @@ dynSparseMatrix::Init_UMFPACK_Sparse(int periods, int y_kmin, int y_kmax, int Si
     }
   else
     {
-      jacob_exo = NULL;
+      jacob_exo = nullptr;
     }
 #ifdef DEBUG
   int local_index;
@@ -1891,7 +1891,7 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair
   line_done = static_cast<bool *>(mxMalloc(Size*periods*sizeof(bool)));
   test_mxMalloc(line_done, __LINE__, __FILE__, __func__, Size*periods*sizeof(bool));
   mem_mngr.init_CHUNK_BLCK_SIZE(u_count);
-  g_save_op = NULL;
+  g_save_op = nullptr;
   g_nop_all = 0;
   i = (periods+y_kmax+1)*Size*sizeof(NonZeroElem *);
   FNZE_R = static_cast<NonZeroElem **>(mxMalloc(i));
@@ -1915,10 +1915,10 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair
     }
   for (int i = 0; i < (periods+y_kmax+1)*Size; i++)
     {
-      FNZE_C[i] = NULL;
-      FNZE_R[i] = NULL;
-      temp_NZE_C[i] = NULL;
-      temp_NZE_R[i] = NULL;
+      FNZE_C[i] = nullptr;
+      FNZE_R[i] = nullptr;
+      temp_NZE_C[i] = nullptr;
+      temp_NZE_R[i] = nullptr;
       NbNZRow[i] = 0;
       NbNZCol[i] = 0;
     }
@@ -1948,19 +1948,19 @@ dynSparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair
                   NbNZRow[eq]++;
                   NbNZCol[var]++;
                   first = mem_mngr.mxMalloc_NZE();
-                  first->NZE_C_N = NULL;
-                  first->NZE_R_N = NULL;
+                  first->NZE_C_N = nullptr;
+                  first->NZE_R_N = nullptr;
                   first->u_index = it4->second+u_count_init*t;
                   first->r_index = eq;
                   first->c_index = var;
                   first->lag_index = lag;
-                  if (FNZE_R[eq] == NULL)
+                  if (FNZE_R[eq] == nullptr)
                     FNZE_R[eq] = first;
-                  if (FNZE_C[var] == NULL)
+                  if (FNZE_C[var] == nullptr)
                     FNZE_C[var] = first;
-                  if (temp_NZE_R[eq] != NULL)
+                  if (temp_NZE_R[eq] != nullptr)
                     temp_NZE_R[eq]->NZE_R_N = first;
-                  if (temp_NZE_C[var] != NULL)
+                  if (temp_NZE_C[var] != nullptr)
                     temp_NZE_C[var]->NZE_C_N = first;
                   temp_NZE_R[eq] = first;
                   temp_NZE_C[var] = first;
@@ -3132,11 +3132,11 @@ dynSparseMatrix::Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned in
     B2_j[++B2_var] = B2_nze;
   while (A3_var < Size)
     A3_j[++A3_var] = A3_nze;
-  mxArray *d1 = NULL;
+  mxArray *d1 = nullptr;
   vector<pair<mxArray *, mxArray *>> triangular_form;
   double sumc = 0, C_sumc = 1000;
-  mxArray *B1_inv = NULL;
-  mxArray *B1_inv_t = NULL;
+  mxArray *B1_inv = nullptr;
+  mxArray *B1_inv_t = nullptr;
   for (int t = 1; t <= periods; t++)
     {
       if (abs(sumc / C_sumc -1) > 1e-10*res1)
@@ -3562,15 +3562,15 @@ dynSparseMatrix::Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double s
 
   status = umfpack_dl_symbolic(n, n, Ap, Ai, Ax, &Symbolic, Control, Info);
   if (status != UMFPACK_OK)
-    umfpack_dl_report_info((double *) NULL, Info);
+    umfpack_dl_report_info(nullptr, Info);
 
   status = umfpack_dl_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info);
   if (status != UMFPACK_OK)
-    umfpack_dl_report_info((double *) NULL, Info);
+    umfpack_dl_report_info(nullptr, Info);
 
   status = umfpack_dl_solve(sys, Ap, Ai, Ax, res, B, Numeric, Control, Info);
   if (status != UMFPACK_OK)
-    umfpack_dl_report_info((double *) NULL, Info);
+    umfpack_dl_report_info(nullptr, Info);
   //double *res = mxGetPr(z);
   if (is_two_boundaries)
     for (int i = 0; i < n; i++)
@@ -4921,9 +4921,9 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
      precond = 1  => Incomplet LU decomposition*/
   size_t n = mxGetM(A_m);
   mxArray *L1, *U1, *Diag;
-  L1 = NULL;
-  U1 = NULL;
-  Diag = NULL;
+  L1 = nullptr;
+  U1 = nullptr;
+  Diag = nullptr;
 
   mxArray *rhs0[4];
   if (preconditioner == 0)
@@ -4978,7 +4978,7 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
     }
   double flags = 2;
   mxArray *z;
-  z = NULL;
+  z = nullptr;
   if (steady_state) /*Octave BicStab algorihtm involves a 0 division in case of a preconditionner equal to the LU decomposition of A matrix*/
     {
       mxArray *res = mult_SAT_B(Sparse_transpose(A_m), x0_m);
@@ -5469,7 +5469,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
 {
   /*Triangularisation at each period of a block using a simple gaussian Elimination*/
   t_save_op_s *save_op_s;
-  int *save_op = NULL, *save_opa = NULL, *save_opaa = NULL;
+  int *save_op = nullptr, *save_opa = nullptr, *save_opaa = nullptr;
   long int nop = 0, nopa = 0;
   bool record = false;
   double *piv_v;
@@ -5983,17 +5983,17 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
                   if (save_opaa)
                     {
                       mxFree(save_opaa);
-                      save_opaa = NULL;
+                      save_opaa = nullptr;
                     }
                   if (save_opa)
                     {
                       mxFree(save_opa);
-                      save_opa = NULL;
+                      save_opa = nullptr;
                     }
                   if (save_op)
                     {
                       mxFree(save_op);
-                      save_op = NULL;
+                      save_op = nullptr;
                     }
                 }
               else if (save_opa && save_opaa)
@@ -6011,7 +6011,7 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
                   if (save_opaa)
                     {
                       mxFree(save_opaa);
-                      save_opaa = NULL;
+                      save_opaa = nullptr;
                     }
                   save_opaa = save_opa;
                 }
@@ -6027,12 +6027,12 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
                   if (save_opa)
                     {
                       mxFree(save_opa);
-                      save_opa = NULL;
+                      save_opa = nullptr;
                     }
                   if (save_opaa)
                     {
                       mxFree(save_opaa);
-                      save_opaa = NULL;
+                      save_opaa = nullptr;
                     }
                 }
             }
@@ -6271,9 +6271,9 @@ bool
 dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, int y_kmax, int size, bool cvg)
 {
   //int i;
-  mxArray *b_m = NULL, *A_m = NULL, *x0_m = NULL;
-  SuiteSparse_long *Ap = NULL, *Ai = NULL;
-  double *Ax = NULL, *b = NULL;
+  mxArray *b_m = nullptr, *A_m = nullptr, *x0_m = nullptr;
+  SuiteSparse_long *Ap = nullptr, *Ai = nullptr;
+  double *Ax = nullptr, *b = nullptr;
   int preconditioner = 1;
 
   try_at_iteration = 0;
@@ -6594,9 +6594,9 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
   u_count_alloc_save = u_count_alloc;
   clock_t t1 = clock();
   nop1 = 0;
-  mxArray *b_m = NULL, *A_m = NULL, *x0_m = NULL;
-  double *Ax = NULL, *b;
-  SuiteSparse_long *Ap = NULL, *Ai = NULL;
+  mxArray *b_m = nullptr, *A_m = nullptr, *x0_m = nullptr;
+  double *Ax = nullptr, *b;
+  SuiteSparse_long *Ap = nullptr, *Ai = nullptr;
 
   if (iter > 0)
     {
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index 1045c81f1..145adb551 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -324,7 +324,7 @@ Get_Arguments_and_global_variables(int nrhs,
                 {
                   *extended_path = true;
                   if ((i+1) >= nrhs)
-                    *ep_struct = NULL;
+                    *ep_struct = nullptr;
                   else
                     {
                       *ep_struct = mxDuplicateArray(prhs[i + 1]);
@@ -369,7 +369,7 @@ Get_Arguments_and_global_variables(int nrhs,
         }
     }
   *M_ = mexGetVariable("global", "M_");
-  if (*M_ == NULL)
+  if (*M_ == nullptr)
     {
       ostringstream tmp;
       tmp << " in main, global variable not found: M_\n";
@@ -377,14 +377,14 @@ Get_Arguments_and_global_variables(int nrhs,
     }
   /* Gets variables and parameters from global workspace of Matlab */
   *oo_ = mexGetVariable("global", "oo_");
-  if (*oo_ == NULL)
+  if (*oo_ == nullptr)
     {
       ostringstream tmp;
       tmp << " in main, global variable not found: oo_\n";
       throw FatalExceptionHandling(tmp.str());
     }
   *options_ = mexGetVariable("global", "options_");
-  if (*options_ == NULL)
+  if (*options_ == nullptr)
     {
       ostringstream tmp;
       tmp << " in main, global variable not found: options_\n";
@@ -404,13 +404,13 @@ main(int nrhs, const char *prhs[])
   mxArray *M_, *oo_, *options_;
   mxArray *GlobalTemporaryTerms;
 #ifndef DEBUG_EX
-  mxArray *block_structur = NULL;
+  mxArray *block_structur = nullptr;
 #else
   int nlhs = 0;
   char *plhs[1];
   load_global((char *) prhs[1]);
 #endif
-  mxArray *pfplan_struct = NULL;
+  mxArray *pfplan_struct = nullptr;
   ErrorMsg error_msg;
   size_t i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
   size_t steady_row_y, steady_col_y;
@@ -420,12 +420,12 @@ main(int nrhs, const char *prhs[])
   bool steady_state = false;
   bool evaluate = false;
   int block = -1;
-  double *params = NULL;
-  double *yd = NULL, *xd = NULL;
+  double *params = nullptr;
+  double *yd = nullptr, *xd = nullptr;
   int count_array_argument = 0;
   bool global_temporary_terms = false;
   bool print = false, print_error = true, print_it = false;
-  double *steady_yd = NULL, *steady_xd = NULL;
+  double *steady_yd = nullptr, *steady_xd = nullptr;
   string plan, pfplan;
   bool extended_path;
   mxArray *extended_path_struct;
@@ -476,13 +476,13 @@ main(int nrhs, const char *prhs[])
 
   if (extended_path)
     {
-      if (extended_path_struct == NULL)
+      if (extended_path_struct == nullptr)
         {
           string tmp = "The 'extended_path' option must be followed by the extended_path descriptor";
           mexErrMsgTxt(tmp.c_str());
         }
       mxArray *date_str = mxGetField(extended_path_struct, 0, "date_str");
-      if (date_str == NULL)
+      if (date_str == nullptr)
         {
           string tmp = "date_str";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
@@ -491,28 +491,28 @@ main(int nrhs, const char *prhs[])
       int nb_periods = mxGetM(date_str) * mxGetN(date_str);
 
       mxArray *constrained_vars_ = mxGetField(extended_path_struct, 0, "constrained_vars_");
-      if (constrained_vars_ == NULL)
+      if (constrained_vars_ == nullptr)
         {
           string tmp = "constrained_vars_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
           mexErrMsgTxt(tmp.c_str());
         }
       mxArray *constrained_paths_ = mxGetField(extended_path_struct, 0, "constrained_paths_");
-      if (constrained_paths_ == NULL)
+      if (constrained_paths_ == nullptr)
         {
           string tmp = "constrained_paths_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
           mexErrMsgTxt(tmp.c_str());
         }
       mxArray *constrained_int_date_ = mxGetField(extended_path_struct, 0, "constrained_int_date_");
-      if (constrained_int_date_ == NULL)
+      if (constrained_int_date_ == nullptr)
         {
           string tmp = "constrained_int_date_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
           mexErrMsgTxt(tmp.c_str());
         }
       mxArray *constrained_perfect_foresight_ = mxGetField(extended_path_struct, 0, "constrained_perfect_foresight_");
-      if (constrained_perfect_foresight_ == NULL)
+      if (constrained_perfect_foresight_ == nullptr)
         {
           string tmp = "constrained_perfect_foresight_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
@@ -520,28 +520,28 @@ main(int nrhs, const char *prhs[])
         }
 
       mxArray *shock_var_ = mxGetField(extended_path_struct, 0, "shock_vars_");
-      if (shock_var_ == NULL)
+      if (shock_var_ == nullptr)
         {
           string tmp = "shock_vars_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
           mexErrMsgTxt(tmp.c_str());
         }
       mxArray *shock_paths_ = mxGetField(extended_path_struct, 0, "shock_paths_");
-      if (shock_paths_ == NULL)
+      if (shock_paths_ == nullptr)
         {
           string tmp = "shock_paths_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
           mexErrMsgTxt(tmp.c_str());
         }
       mxArray *shock_int_date_ = mxGetField(extended_path_struct, 0, "shock_int_date_");
-      if (shock_int_date_ == NULL)
+      if (shock_int_date_ == nullptr)
         {
           string tmp = "shock_int_date_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
           mexErrMsgTxt(tmp.c_str());
         }
       mxArray *shock_str_date_ = mxGetField(extended_path_struct, 0, "shock_str_date_");
-      if (shock_str_date_ == NULL)
+      if (shock_str_date_ == nullptr)
         {
           string tmp = "shock_str_date_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
@@ -550,8 +550,8 @@ main(int nrhs, const char *prhs[])
       int nb_constrained = mxGetM(constrained_vars_) * mxGetN(constrained_vars_);
       int nb_controlled = 0;
       mxArray *options_cond_fcst_ = mxGetField(extended_path_struct, 0, "options_cond_fcst_");
-      mxArray *controlled_varexo = NULL;
-      if (options_cond_fcst_ != NULL)
+      mxArray *controlled_varexo = nullptr;
+      if (options_cond_fcst_ != nullptr)
         {
           controlled_varexo = mxGetField(options_cond_fcst_, 0, "controlled_varexo");
           nb_controlled = mxGetM(controlled_varexo) * mxGetN(controlled_varexo);
@@ -560,8 +560,8 @@ main(int nrhs, const char *prhs[])
               mexErrMsgTxt("The number of exogenized variables and the number of exogenous controlled variables should be equal.");
             }
         }
-      double *controlled_varexo_value = NULL;
-      if (controlled_varexo != NULL)
+      double *controlled_varexo_value = nullptr;
+      if (controlled_varexo != nullptr)
         controlled_varexo_value = mxGetPr(controlled_varexo);
       double *constrained_var_value = mxGetPr(constrained_vars_);
       sconditional_extended_path.resize(nb_constrained);
@@ -680,7 +680,7 @@ main(int nrhs, const char *prhs[])
   if (plan.length() > 0)
     {
       mxArray *plan_struct = mexGetVariable("base", plan.c_str());
-      if (plan_struct == NULL)
+      if (plan_struct == nullptr)
         {
           string tmp = plan;
           tmp.insert(0, "Can't find the plan: ");