diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 793afa56164ad0fe4f39aa2f7f61b14b0a318617..ee1546af50d5fe280641e199cfcc0a681cf5532f 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -58,7 +58,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
 % SPECIAL REQUIREMENTS
 %   None
 
-% Copyright (C) 2006-2018 Dynare Team
+% Copyright (C) 2006-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -185,8 +185,7 @@ elseif options_.lik_init == 3           % Diffuse Kalman filter
     end
     [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium);
 elseif options_.lik_init == 4           % Start from the solution of the Riccati equation.
-    [err, Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,vobs)),H);
-    mexErrCheck('kalman_steady_state',err);
+    Pstar = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,vobs)),H);
     Pinf  = [];
     if kalman_algo~=2
         kalman_algo = 1;
diff --git a/matlab/block_bytecode_mfs_steadystate.m b/matlab/block_bytecode_mfs_steadystate.m
index a34b15c197196a4d31dd342812c36d76856b14e8..ba0bef8617b4d48e2aefd7b85bfaff6c36480588 100644
--- a/matlab/block_bytecode_mfs_steadystate.m
+++ b/matlab/block_bytecode_mfs_steadystate.m
@@ -2,7 +2,7 @@ function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, M)
 % Wrapper around the *_static.m file, for use with dynare_solve,
 % when block_mfs option is given to steady.
 
-% Copyright (C) 2009-2012 Dynare Team
+% Copyright (C) 2009-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -21,4 +21,4 @@ function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, M)
 
 indx = M.block_structure_stat.block(b).variable;
 y_all(indx) = y;
-[chk, r, g1] = bytecode( y_all, exo, params, y_all, 1, y_all, 'evaluate', 'static', ['block = ' int2str(b) ]);
+[r, g1] = bytecode( y_all, exo, params, y_all, 1, y_all, 'evaluate', 'static', ['block = ' int2str(b) ]);
diff --git a/matlab/bytecode_steadystate.m b/matlab/bytecode_steadystate.m
index 4c3e2f46f9e4883d906171c038add32e5e15ada7..dff35373f4875072b510294ab2d7c97f84e84180 100644
--- a/matlab/bytecode_steadystate.m
+++ b/matlab/bytecode_steadystate.m
@@ -2,7 +2,7 @@ function [r, g1] = bytecode_steadystate(y, exo, params)
 % Wrapper around the *_static.m file, for use with dynare_solve,
 % when block_mfs option is given to steady.
 
-% Copyright (C) 2009-2011 Dynare Team
+% Copyright (C) 2009-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -19,4 +19,4 @@ function [r, g1] = bytecode_steadystate(y, exo, params)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-eval('[chk, r, g1] = bytecode( y, exo, params, y, 1, exo, ''evaluate'', ''static'');');
+eval('[r, g1] = bytecode( y, exo, params, y, 1, exo, ''evaluate'', ''static'');');
diff --git a/matlab/dr_block.m b/matlab/dr_block.m
index 0328c1048c24ec3432562717e137cc332a2eed34..d65248b7745ccd275598803f94d790888d7fe0cb 100644
--- a/matlab/dr_block.m
+++ b/matlab/dr_block.m
@@ -34,7 +34,7 @@ function [dr,info,M_,options_,oo_] = dr_block(dr,task,M_,options_,oo_,varargin)
 %   none.
 %
 
-% Copyright (C) 2010-2017 Dynare Team
+% Copyright (C) 2010-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -71,12 +71,10 @@ else
     Size = 1;
 end
 if (options_.bytecode)
-    [chck, zz, data]= bytecode('dynamic','evaluate', z, zx, M_.params, dr.ys, 1, data);
+    [zz, data]= bytecode('dynamic','evaluate', z, zx, M_.params, dr.ys, 1, data);
 else
     [r, data] = feval([M_.fname '.dynamic'], options_, M_, oo_, z', zx, M_.params, dr.ys, M_.maximum_lag+1, data);
-    chck = 0;
 end
-mexErrCheck('bytecode', chck);
 dr.full_rank = 1;
 dr.eigval = [];
 dr.nd = 0;
@@ -440,7 +438,7 @@ for i = 1:Size
             D = [[aa(row_indx,index_0m) zeros(n_dynamic,n_both) aa(row_indx,index_p)] ; [zeros(n_both, n_pred) eye(n_both) zeros(n_both, n_both + n_fwrd)]];
             E = [-aa(row_indx,[index_m index_0p])  ; [zeros(n_both, n_both + n_pred) eye(n_both, n_both + n_fwrd) ] ];
 
-            [err, ss, tt, w, sdim, data(i).eigval, info1] = mjdgges(E,D,options_.qz_criterium,options_.qz_zero_threshold);
+            [ss, tt, w, sdim, data(i).eigval, info1] = mjdgges(E,D,options_.qz_criterium,options_.qz_zero_threshold);
 
             if (verbose)
                 disp('eigval');
@@ -591,7 +589,7 @@ for i = 1:Size
                 elseif options_.sylvester_fp
                     ghx_other = gensylv_fp(A_, B_, C_, D_, i, options_.sylvester_fixed_point_tol);
                 else
-                    [err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
+                    ghx_other = gensylv(1, A_, B_, C_, -D_);
                 end
                 if options_.aim_solver ~= 1
                     % Necessary when using Sims' routines for QZ
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 40a7203c43bccbadad3245141b5f648dd480bd08..ddcdf33f40db400399a0b0d7e5f85d8589ce6cb7 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -115,7 +115,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 2004-2019 Dynare Team
+% Copyright (C) 2004-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -150,7 +150,9 @@ if DynareOptions.estimation_dll
     [fval,exit_flag,SteadyState,trend_coeff,info,params,H,Q] ...
         = logposterior(xparam1,DynareDataset, DynareOptions,Model, ...
                        EstimatedParameters,BayesInfo,DynareResults);
-    mexErrCheck('logposterior', exit_flag);
+    if exit_flag
+        error("Error encountered in logposterior")
+    end
     Model.params = params;
     if ~isequal(Model.H,0)
         Model.H = H;
@@ -455,12 +457,14 @@ switch DynareOptions.lik_init
     if kalman_algo ~= 2
         kalman_algo = 1;
     end
-    if isequal(H,0)
-        [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,mm,length(Z))));
-    else
-        [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,mm,length(Z))),H);
-    end
-    if err
+    try
+        if isequal(H,0)
+            Pstar = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,mm,length(Z))));
+        else
+            Pstar = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,mm,length(Z))),H);
+        end
+    catch ME
+        disp(ME.message)
         disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
         DynareOptions.lik_init = 1;
         Pstar=lyapunov_solver(T,R,Q,DynareOptions);
@@ -647,8 +651,7 @@ singularity_has_been_detected = false;
 if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
     if no_missing_data_flag
         if DynareOptions.block
-            [err, LIK] = block_kalman_filter(T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
-            mexErrCheck('block_kalman_filter', err);
+            LIK = block_kalman_filter(T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
         elseif DynareOptions.fast_kalman_filter
             if diffuse_periods
                 %kalman_algo==3 requires no diffuse periods (stationary
@@ -677,8 +680,8 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
         end
     else
         if 0 %DynareOptions.block
-            [err, LIK,lik] = block_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,...
-                                                 T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
+            [LIK,lik] = block_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,...
+                                            T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
         else
             [LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
                                                            a, Pstar, ...
diff --git a/matlab/dyn_first_order_solver.m b/matlab/dyn_first_order_solver.m
index 63cb4f9da5a00fd13edcb9935b76f74a3d819b3b..101a33700c03032b7955e417ae8a1f77ab867b02 100644
--- a/matlab/dyn_first_order_solver.m
+++ b/matlab/dyn_first_order_solver.m
@@ -21,7 +21,7 @@ function [dr, info] = dyn_first_order_solver(jacobia, DynareModel, dr, DynareOpt
 %                                     info=5 -> Blanchard and Kahn conditions are not satisfied: indeterminacy due to rank failure,
 %                                     info=7 -> One of the eigenvalues is close to 0/0 (infinity of complex solutions)
 
-% Copyright (C) 2001-2018 Dynare Team
+% Copyright (C) 2001-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -187,8 +187,7 @@ else
     E(row_indx_de_1,index_e1) = -aa(row_indx,index_e);
     E(row_indx_de_2,index_e2) = eye(nboth);
 
-    [err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium, DynareOptions.qz_zero_threshold);
-    mexErrCheck('mjdgges', err);
+    [ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium, DynareOptions.qz_zero_threshold);
 
     if info1
         if info1 == -30
diff --git a/matlab/dyn_ramsey_static.m b/matlab/dyn_ramsey_static.m
index e94deb5044772f6df93fd26f244be6269ed430e0..6ea85e7d6e410b61a5c4c9720e653735c9a5f84c 100644
--- a/matlab/dyn_ramsey_static.m
+++ b/matlab/dyn_ramsey_static.m
@@ -18,7 +18,7 @@ function [steady_state,params,check] = dyn_ramsey_static(ys_init,M,options_,oo)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2003-2018 Dynare Team
+% Copyright (C) 2003-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -156,8 +156,8 @@ Uyy = reshape(Uyy,endo_nbr,endo_nbr);
 % set multipliers and auxiliary variables that
 % depends on multipliers to 0 to compute residuals
 if (options_.bytecode)
-    [chck, res, junk] = bytecode('static',xx,[oo.exo_steady_state oo.exo_det_steady_state], ...
-                                 params, 'evaluate');
+    [res, junk] = bytecode('static',xx,[oo.exo_steady_state oo.exo_det_steady_state], ...
+                           params, 'evaluate');
     fJ = junk.g1;
 else
     [res,fJ] = feval([fname '.static'],xx,[oo.exo_steady_state oo.exo_det_steady_state], ...
@@ -194,8 +194,8 @@ end
 function result = check_static_model(ys,M,options_,oo)
 result = false;
 if (options_.bytecode)
-    [chck, res, ~] = bytecode('static',ys,[oo.exo_steady_state oo.exo_det_steady_state], ...
-                              M.params, 'evaluate');
+    [res, ~] = bytecode('static',ys,[oo.exo_steady_state oo.exo_det_steady_state], ...
+                        M.params, 'evaluate');
 else
     res = feval([M.fname '.static'],ys,[oo.exo_steady_state oo.exo_det_steady_state], ...
                 M.params);
diff --git a/matlab/dyn_risky_steadystate_solver.m b/matlab/dyn_risky_steadystate_solver.m
index d8fe2ede16a647a3edae27e3baaa42ddccc8589a..93c75da4dc4d2b98b92f3f9fafca2668fe0b888e 100644
--- a/matlab/dyn_risky_steadystate_solver.m
+++ b/matlab/dyn_risky_steadystate_solver.m
@@ -341,8 +341,7 @@ if nargout > 1
     nu2 = exo_nbr*(exo_nbr+1)/2;
     nu3 = exo_nbr*(exo_nbr+1)*(exo_nbr+2)/3;
     M_np.NZZDerivatives = [nnz(d1_np); nnz(d2_np); nnz(d3_np)];
-    [err, dynpp_derivs] = k_order_perturbation(dr_np,M_np,options,d1_np,d2_np,d3_np);
-    mexErrCheck('k_order_perturbation', err);
+    dynpp_derivs = k_order_perturbation(dr_np,M_np,options,d1_np,d2_np,d3_np);
     g_0 = dynpp_derivs.g_0;
     g_1 = dynpp_derivs.g_1;
     g_2 = dynpp_derivs.g_2;
diff --git a/matlab/dyn_second_order_solver.m b/matlab/dyn_second_order_solver.m
index 4e3151497c1f0f0641d5c6bc6792cbfd1fd33481..bd5c37a9fa786dd19eac24ef12e396ff3563e658 100644
--- a/matlab/dyn_second_order_solver.m
+++ b/matlab/dyn_second_order_solver.m
@@ -36,7 +36,7 @@ function dr = dyn_second_order_solver(jacobia,hessian_mat,dr,M,threads_BC)
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 2001-2019 Dynare Team
+% Copyright (C) 2001-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -85,29 +85,23 @@ zu = [zeros(length(ic),M.exo_nbr);
       dr.ghx(klead~=0,:)*dr.ghu(ic,:);
       eye(M.exo_nbr);
       zeros(M.exo_det_nbr,M.exo_nbr)];
-[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,threads_BC); %hessian_mat: reordering to DR order
-mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
+rhs = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,threads_BC); %hessian_mat: reordering to DR order
 rhs = -rhs;
-[err, dr.ghxx] = gensylv(2,A,B,C,rhs);
-mexErrCheck('gensylv', err);
+dr.ghxx = gensylv(2,A,B,C,rhs);
 
 
 %% ghxu
 %rhs
-[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,zu,threads_BC); %hessian_mat: reordering to DR order
-mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-[abcOut,err] = A_times_B_kronecker_C(dr.ghxx, dr.ghx(ic,:), dr.ghu(ic,:));
-mexErrCheck('A_times_B_kronecker_C', err);
+rhs = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,zu,threads_BC); %hessian_mat: reordering to DR order
+abcOut = A_times_B_kronecker_C(dr.ghxx, dr.ghx(ic,:), dr.ghu(ic,:));
 rhs = -rhs-B*abcOut;
 %lhs
 dr.ghxu = A\rhs;
 
 %% ghuu
 %rhs
-[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zu,threads_BC); %hessian_mat: reordering to DR order
-mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-[B1, err] = A_times_B_kronecker_C(B*dr.ghxx,dr.ghu(ic,:));
-mexErrCheck('A_times_B_kronecker_C', err);
+rhs = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zu,threads_BC); %hessian_mat: reordering to DR order
+B1 = A_times_B_kronecker_C(B*dr.ghxx,dr.ghu(ic,:));
 rhs = -rhs-B1;
 %lhs
 dr.ghuu = A\rhs;
@@ -120,8 +114,7 @@ LHS = zeros(M.endo_nbr,M.endo_nbr);
 LHS(:,kcurr~=0) = jacobia(:,nonzeros(kcurr));
 RHS = zeros(M.endo_nbr,M.exo_nbr^2);
 E = eye(M.endo_nbr);
-[B1, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(nonzeros(klead),nonzeros(klead))), dr.ghu(klead~=0,:),threads_BC); %hessian_mat:focus only on forward variables and reorder to DR order
-mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
+B1 = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(nonzeros(klead),nonzeros(klead))), dr.ghu(klead~=0,:),threads_BC); %hessian_mat:focus only on forward variables and reorder to DR order
 RHS = RHS + jacobia(:,nonzeros(klead))*dr.ghuu(klead~=0,:)+B1;
 % LHS
 LHS = LHS + jacobia(:,nonzeros(klead))*(E(klead~=0,:)+[O1(klead~=0,:) dr.ghx(klead~=0,:) O2(klead~=0,:)]);
diff --git a/matlab/dynare_solve_block_or_bytecode.m b/matlab/dynare_solve_block_or_bytecode.m
index c995c93c7d7d3d3d337c0ab40c4e62bb2c6dd5c7..f991f54f53e725f248fb0e5244fcbaa2581265b3 100644
--- a/matlab/dynare_solve_block_or_bytecode.m
+++ b/matlab/dynare_solve_block_or_bytecode.m
@@ -1,5 +1,5 @@
 function [x,info] = dynare_solve_block_or_bytecode(y, exo, params, options, M)
-% Copyright (C) 2010-2017 Dynare Team
+% Copyright (C) 2010-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -51,9 +51,10 @@ if options.block && ~options.bytecode
     end
 elseif options.bytecode
     if options.solve_algo > 4
-        [check, x] = bytecode('static', x, exo, params);
-        %        mexErrCheck('bytecode', check);
-        if check
+        try
+            x = bytecode('static', x, exo, params);
+        catch ME
+            disp(ME.message);
             info = 1;
             return
         end
@@ -73,10 +74,12 @@ elseif options.bytecode
                 end
                 x(M.block_structure_stat.block(b).variable) = y;
             else
-                [chk, nulldev, nulldev1, x] = bytecode( x, exo, params, ...
-                                                        x, 1, x, 'evaluate', 'static', ...
-                                                        ['block = ' int2str(b)]);
-                if chk
+                try
+                    [nulldev, nulldev1, x] = bytecode(x, exo, params, ...
+                                                      x, 1, x, 'evaluate', 'static', ...
+                                                      ['block = ' int2str(b)]);
+                catch ME
+                    disp(ME.message);
                     info = 1;
                     return
                 end
diff --git a/matlab/ep/extended_path_core.m b/matlab/ep/extended_path_core.m
index 1768c0bbd60731f0813270540c766c37c6497d4b..eb8b2454dd7d0a220c2d35443994491be8147654 100644
--- a/matlab/ep/extended_path_core.m
+++ b/matlab/ep/extended_path_core.m
@@ -4,7 +4,7 @@ function [y, info_convergence, endogenousvariablespaths] = extended_path_core(pe
                                                   debug,bytecode_flag,order,M,pfm,algo,solve_algo,stack_solve_algo,...
                                                   olmmcp,options,oo,initialguess)
 
-% Copyright (C) 2016-2019 Dynare Team
+% Copyright (C) 2016-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -41,7 +41,13 @@ if debug
 end
 
 if bytecode_flag && ~ep.stochastic.order
-    [flag, tmp] = bytecode('dynamic', endo_simul, exo_simul, M_.params, endo_simul, periods);
+    try
+        tmp = bytecode('dynamic', endo_simul, exo_simul, M_.params, endo_simul, periods);
+        flag = false;
+    catch ME
+        disp(ME.message);
+        flag = true;
+    end
 else
     flag = true;
 end
@@ -86,4 +92,4 @@ else
     y = NaN(size(endo_nbr,1));
 end
 
-endogenousvariablespaths = tmp.endo_simul;
\ No newline at end of file
+endogenousvariablespaths = tmp.endo_simul;
diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m
index 198d65ccda5dc246519fdf830f3fe7b947603172..63ec0ab6a7c380c7a0829a54fd03e8d755e00353 100644
--- a/matlab/evaluate_planner_objective.m
+++ b/matlab/evaluate_planner_objective.m
@@ -11,7 +11,7 @@ function planner_objective_value = evaluate_planner_objective(M,options,oo)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2007-2018 Dynare Team
+% Copyright (C) 2007-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -57,21 +57,16 @@ ys = oo.dr.ys;
 %second order terms
 Uyy = full(Uyy);
 
-[Uyygygy, err] = A_times_B_kronecker_C(Uyy,gy,gy);
-mexErrCheck('A_times_B_kronecker_C', err);
-[Uyygugu, err] = A_times_B_kronecker_C(Uyy,gu,gu);
-mexErrCheck('A_times_B_kronecker_C', err);
-[Uyygygu, err] = A_times_B_kronecker_C(Uyy,gy,gu);
-mexErrCheck('A_times_B_kronecker_C', err);
+Uyygygy = A_times_B_kronecker_C(Uyy,gy,gy);
+Uyygugu = A_times_B_kronecker_C(Uyy,gu,gu);
+Uyygygu = A_times_B_kronecker_C(Uyy,gy,gu);
 
 Wbar =U/(1-beta); %steady state welfare
 Wy = Uy*gy/(eye(nspred)-beta*Gy);
 Wu = Uy*gu+beta*Wy*Gu;
 Wyy = Uyygygy/(eye(nspred*nspred)-beta*kron(Gy,Gy));
-[Wyygugu, err] = A_times_B_kronecker_C(Wyy,Gu,Gu);
-mexErrCheck('A_times_B_kronecker_C', err);
-[Wyygygu,err] = A_times_B_kronecker_C(Wyy,Gy,Gu);
-mexErrCheck('A_times_B_kronecker_C', err);
+Wyygugu = A_times_B_kronecker_C(Wyy,Gu,Gu);
+Wyygygu = A_times_B_kronecker_C(Wyy,Gy,Gu);
 Wuu = Uyygugu+beta*Wyygugu;
 Wyu = Uyygygu+beta*Wyygygu;
 Wss = beta*Wuu*M.Sigma_e(:)/(1-beta); % at period 0, we are in steady state, so the deviation term only starts in period 1, thus the beta in front
@@ -93,20 +88,15 @@ end
 yhat1 = yhat1(dr.order_var(nstatic+(1:nspred)),1)-dr.ys(dr.order_var(nstatic+(1:nspred)));
 u = oo.exo_simul(1,:)';
 
-[Wyyyhatyhat1, err] = A_times_B_kronecker_C(Wyy,yhat1,yhat1);
-mexErrCheck('A_times_B_kronecker_C', err);
-[Wuuuu, err] = A_times_B_kronecker_C(Wuu,u,u);
-mexErrCheck('A_times_B_kronecker_C', err);
-[Wyuyhatu1, err] = A_times_B_kronecker_C(Wyu,yhat1,u);
-mexErrCheck('A_times_B_kronecker_C', err);
+Wyyyhatyhat1 = A_times_B_kronecker_C(Wyy,yhat1,yhat1);
+Wuuuu = A_times_B_kronecker_C(Wuu,u,u);
+Wyuyhatu1 = A_times_B_kronecker_C(Wyu,yhat1,u);
 planner_objective_value(1) = Wbar+Wy*yhat1+Wu*u+Wyuyhatu1 ...
     + 0.5*(Wyyyhatyhat1 + Wuuuu+Wss);
 if options.ramsey_policy
     yhat2 = yhat2(dr.order_var(nstatic+(1:nspred)),1)-dr.ys(dr.order_var(nstatic+(1:nspred)));
-    [Wyyyhatyhat2, err] = A_times_B_kronecker_C(Wyy,yhat2,yhat2);
-    mexErrCheck('A_times_B_kronecker_C', err);
-    [Wyuyhatu2, err] = A_times_B_kronecker_C(Wyu,yhat2,u);
-    mexErrCheck('A_times_B_kronecker_C', err);
+    Wyyyhatyhat2 = A_times_B_kronecker_C(Wyy,yhat2,yhat2);
+    Wyuyhatu2 = A_times_B_kronecker_C(Wyu,yhat2,u);
     planner_objective_value(2) = Wbar+Wy*yhat2+Wu*u+Wyuyhatu2 ...
         + 0.5*(Wyyyhatyhat2 + Wuuuu+Wss);
 end
diff --git a/matlab/evaluate_static_model.m b/matlab/evaluate_static_model.m
index d15005aa208eaa473b4629e712903c5d558919bb..fbcfa235cd30edc8b3b68cc19209b96ae93d6030 100644
--- a/matlab/evaluate_static_model.m
+++ b/matlab/evaluate_static_model.m
@@ -20,7 +20,7 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2001-2017 Dynare Team
+% Copyright (C) 2001-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -39,9 +39,8 @@ function [residuals,check1,jacob] = evaluate_static_model(ys,exo_ss,params,M,opt
 
 check1 = 0;
 if options.bytecode
-    [check1, residuals] = bytecode('evaluate','static',ys,...
-                                   exo_ss, params, ys, 1);
-    mexErrCheck('bytecode', check1);
+    residuals = bytecode('evaluate','static',ys,...
+                         exo_ss, params, ys, 1);
 else
     fh_static = str2func([M.fname '.static']);
     if options.block
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 12507a2b4686e964a6707b1668e3316553015d68..8fc5826c875e95af34205bb8fe0ccc96436a0a1e 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -22,7 +22,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2001-2018 Dynare Team
+% Copyright (C) 2001-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -307,8 +307,7 @@ if M.static_and_dynamic_models_differ
     z = repmat(ys,1,M.maximum_lead + M.maximum_lag + 1);
     zx = repmat([exo_ss'], M.maximum_lead + M.maximum_lag + 1, 1);
     if options.bytecode
-        [chck, r, ~]= bytecode('dynamic','evaluate', z, zx, params, ys, 1);
-        mexErrCheck('bytecode', chck);
+        [r, ~]= bytecode('dynamic','evaluate', z, zx, params, ys, 1);
     elseif options.block
         [r, oo.dr] = feval([M.fname '.dynamic'], z', zx, params, ys, M.maximum_lag+1, oo.dr);
     else
diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/get_perturbation_params_derivs.m
index 64c5a2e3469f3a55d9428784047dc0d2d44dd385..52b5f9cecf450b46f73751235fd09c15a4d1cca9 100644
--- a/matlab/get_perturbation_params_derivs.m
+++ b/matlab/get_perturbation_params_derivs.m
@@ -97,7 +97,7 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i
 %   * sylvester3a
 %   * get_perturbation_params_derivs_numerical_objective
 % =========================================================================
-% Copyright (C) 2019 Dynare Team
+% Copyright (C) 2019-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -821,9 +821,9 @@ if order > 1
                        dghx(kcurr~=0,:,jp);
                        dghx(klead~=0,:,jp)*oo.dr.ghx(idx_states,:) + oo.dr.ghx(klead~=0,:)*dghx(idx_states,:,jp);
                        zeros(M.exo_nbr,M.nspred)];
-        [dRHS_1, err] = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0,kyy0)+(jp-1)*yy0ex0_nbr^2),zx(1:yy0_nbr,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dRHS_2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0,kyy0)),dzx(1:yy0_nbr,:,jp),zx(1:yy0_nbr,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dRHS_3, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0,kyy0)),zx(1:yy0_nbr,:),dzx(1:yy0_nbr,:,jp),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
+        dRHS_1 = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0,kyy0)+(jp-1)*yy0ex0_nbr^2),zx(1:yy0_nbr,:),threads_BC);
+        dRHS_2 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0,kyy0)),dzx(1:yy0_nbr,:,jp),zx(1:yy0_nbr,:),threads_BC);
+        dRHS_3 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0,kyy0)),zx(1:yy0_nbr,:),dzx(1:yy0_nbr,:,jp),threads_BC);
         dG_x_x(:,:,jp) = kron(dghx(idx_states,:,jp),oo.dr.ghx(idx_states,:)) + kron(oo.dr.ghx(idx_states,:),dghx(idx_states,:,jp));
         dRHSghxx(:,:,jp) = -( (dRHS_1+dRHS_2+dRHS_3) + dA(:,:,jp)*oo.dr.ghxx + dB(:,:,jp)*oo.dr.ghxx*G_x_x + B*oo.dr.ghxx*dG_x_x(:,:,jp) );
     end
@@ -838,9 +838,9 @@ if order > 1
           oo.dr.ghu(kcurr~=0,:);
           oo.dr.ghx(klead~=0,:)*oo.dr.ghu(idx_states,:);
           eye(M.exo_nbr)];
-    [abcOutxu,err] = A_times_B_kronecker_C(oo.dr.ghxx,oo.dr.ghx(idx_states,:),oo.dr.ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err); %auxiliary expressions for dghxu
-    [abcOutuu, err] = A_times_B_kronecker_C(oo.dr.ghxx,oo.dr.ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err); %auxiliary expressions for dghuu
-    [RHSs2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), oo.dr.ghu(klead~=0,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
+    abcOutxu = A_times_B_kronecker_C(oo.dr.ghxx,oo.dr.ghx(idx_states,:),oo.dr.ghu(idx_states,:)); %auxiliary expressions for dghxu
+    abcOutuu = A_times_B_kronecker_C(oo.dr.ghxx,oo.dr.ghu(idx_states,:)); %auxiliary expressions for dghuu
+    RHSs2 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), oo.dr.ghu(klead~=0,:),threads_BC);
     RHSs2 = RHSs2 + g1(:,nonzeros(klead))*oo.dr.ghuu(klead~=0,:);
     dzu = zeros(size(zu,1),size(zu,2),modparam_nbr);
     dghxu = zeros(M.endo_nbr,M.nspred*M.exo_nbr,modparam_nbr);
@@ -852,25 +852,25 @@ if order > 1
                        dghx(klead~=0,:,jp)*oo.dr.ghu(idx_states,:) + oo.dr.ghx(klead~=0,:)*dghu(idx_states,:,jp);
                        zeros(M.exo_nbr,M.exo_nbr)];
         %compute dghxu
-        [dRHS_1, err] = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2),zx,zu,threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dRHS_2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),dzx(:,:,jp),zu,threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dRHS_3, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),zx,dzu(:,:,jp),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dabcOut_1,err] = A_times_B_kronecker_C(dghxx(:,:,jp),oo.dr.ghx(idx_states,:),oo.dr.ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err);
-        [dabcOut_2,err] = A_times_B_kronecker_C(oo.dr.ghxx,dghx(idx_states,:,jp),oo.dr.ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err);
-        [dabcOut_3,err] = A_times_B_kronecker_C(oo.dr.ghxx,oo.dr.ghx(idx_states,:),dghu(idx_states,:,jp)); mexErrCheck('A_times_B_kronecker_C', err);
+        dRHS_1 = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2),zx,zu,threads_BC);
+        dRHS_2 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),dzx(:,:,jp),zu,threads_BC);
+        dRHS_3 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),zx,dzu(:,:,jp),threads_BC);
+        dabcOut_1 = A_times_B_kronecker_C(dghxx(:,:,jp),oo.dr.ghx(idx_states,:),oo.dr.ghu(idx_states,:));
+        dabcOut_2 = A_times_B_kronecker_C(oo.dr.ghxx,dghx(idx_states,:,jp),oo.dr.ghu(idx_states,:));
+        dabcOut_3 = A_times_B_kronecker_C(oo.dr.ghxx,oo.dr.ghx(idx_states,:),dghu(idx_states,:,jp));
         dghxu(:,:,jp) = -invA * ( dA(:,:,jp)*oo.dr.ghxu + (dRHS_1+dRHS_2+dRHS_3) + dB(:,:,jp)*abcOutxu + B*(dabcOut_1+dabcOut_2+dabcOut_3) );
         %compute dghuu
-        [dRHS_1, err] = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2),zu,threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dRHS_2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),dzu(:,:,jp),zu,threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dRHS_3, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),zu,dzu(:,:,jp),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dabcOut_1, err] = A_times_B_kronecker_C(dghxx(:,:,jp),oo.dr.ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err);
-        [dabcOut_2, err] = A_times_B_kronecker_C(oo.dr.ghxx,dghu(idx_states,:,jp),oo.dr.ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err);
-        [dabcOut_3, err] = A_times_B_kronecker_C(oo.dr.ghxx,oo.dr.ghu(idx_states,:),dghu(idx_states,:,jp)); mexErrCheck('A_times_B_kronecker_C', err);
+        dRHS_1 = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2),zu,threads_BC);
+        dRHS_2 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),dzu(:,:,jp),zu,threads_BC);
+        dRHS_3 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),zu,dzu(:,:,jp),threads_BC);
+        dabcOut_1 = A_times_B_kronecker_C(dghxx(:,:,jp),oo.dr.ghu(idx_states,:));
+        dabcOut_2 = A_times_B_kronecker_C(oo.dr.ghxx,dghu(idx_states,:,jp),oo.dr.ghu(idx_states,:));
+        dabcOut_3 = A_times_B_kronecker_C(oo.dr.ghxx,oo.dr.ghu(idx_states,:),dghu(idx_states,:,jp));
         dghuu(:,:,jp) = -invA * ( dA(:,:,jp)*oo.dr.ghuu + (dRHS_1+dRHS_2+dRHS_3) + dB(:,:,jp)*abcOutuu + B*(dabcOut_1+dabcOut_2+dabcOut_3) );
         %compute dghs2
-        [dRHS_1, err] = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))+(jp-1)*yy0ex0_nbr^2), oo.dr.ghu(klead~=0,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dRHS_2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), dghu(klead~=0,:,jp), oo.dr.ghu(klead~=0,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
-        [dRHS_3, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), oo.dr.ghu(klead~=0,:), dghu(klead~=0,:,jp),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
+        dRHS_1 = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))+(jp-1)*yy0ex0_nbr^2), oo.dr.ghu(klead~=0,:),threads_BC);
+        dRHS_2 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), dghu(klead~=0,:,jp), oo.dr.ghu(klead~=0,:),threads_BC);
+        dRHS_3 = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), oo.dr.ghu(klead~=0,:), dghu(klead~=0,:,jp),threads_BC);
         dghs2(:,stderrparam_nbr+corrparam_nbr+jp) = -invS * ( dS(:,:,jp)*oo.dr.ghs2 + ( (dRHS_1+dRHS_2+dRHS_3) + dg1(:,nonzeros(klead),jp)*oo.dr.ghuu(klead~=0,:) + g1(:,nonzeros(klead))*dghuu(klead~=0,:,jp) )*M.Sigma_e(:) );
     end
     %add contributions by dSigma_e to dghs2
diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m
index 7627970cec4ecde3d05bdcf60270902e01ace885..beaa621b68212ec40d60aa1ed6ac72497d062dc9 100644
--- a/matlab/k_order_pert.m
+++ b/matlab/k_order_pert.m
@@ -1,7 +1,7 @@
 function [dr,info] = k_order_pert(dr,M,options)
 % Compute decision rules using the k-order DLL from Dynare++
 
-% Copyright (C) 2009-2018 Dynare Team
+% Copyright (C) 2009-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -30,8 +30,9 @@ if M.maximum_endo_lead == 0 && order>1
            'backward models'])
 end
 
-[err, dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,M,options);
-if err
+try
+    [dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,M,options);
+catch
     info(1)=9;
     return
 end
diff --git a/matlab/mex/k_order_perturbation.m b/matlab/mex/k_order_perturbation.m
index 2f4b1eec96723d045e5b688d6571415b91a2bc81..15e34a7ce86f73925205ff705e8111fc5333062f 100644
--- a/matlab/mex/k_order_perturbation.m
+++ b/matlab/mex/k_order_perturbation.m
@@ -1,4 +1,4 @@
-% [err, dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,DynareModel,DynareOptions,g1,g2,g3)
+% [dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,DynareModel,DynareOptions,g1,g2,g3)
 % computes a k_order_petrubation solution for k=1,2,3
 %
 % INPUTS
@@ -10,7 +10,6 @@
 % g3:            matrix   Third order derivatives of the model (optional)
 %
 % OUTPUTS
-% err:           double   error code
 % dynpp_derivs   struct   Derivatives of the decision rule in Dynare++ format.
 %                         The tensors are folded and the Taylor coefficients (1/n!) are
 %                         included.
@@ -29,7 +28,7 @@
 % dynare/mex/sources/k_order_perturbation.cc and it uses code provided by
 % dynare++
 
-% Copyright (C) 2013-2017 Dynare Team
+% Copyright (C) 2013-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
diff --git a/matlab/mexErrCheck.m b/matlab/mexErrCheck.m
deleted file mode 100644
index 708d4e5057205fed732648286d4fa58b7cac892d..0000000000000000000000000000000000000000
--- a/matlab/mexErrCheck.m
+++ /dev/null
@@ -1,42 +0,0 @@
-function mexErrCheck(mexFunctionName, err)
-% function mexErrCheck(mexFunctionName, err)
-% this function halts processing if err is equal to 1.
-%
-% INPUTS
-%   mexFunctionName [char]    Name of the mexFunction
-%   err             [double]  error code returned from mexFunction
-%
-% OUTPUTS
-%   none.
-%
-% ALGORITHM
-%   ...
-%
-% SPECIAL REQUIREMENTS
-%   none.
-%
-
-% Copyright (C) 2010 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
-
-if ~ischar(mexFunctionName) || ~isscalar(err)
-    error('The first argument must be a char and the second a scalar');
-end
-
-if err
-    error(['Error encountered in: ' mexFunctionName '.']);
-end
\ No newline at end of file
diff --git a/matlab/missing/mex/gensylv/gensylv.m b/matlab/missing/mex/gensylv/gensylv.m
index 2920aca4d2bd8d3f6f6cfb0842618cf669e75673..462e7b1ea0f0eaecaa1dc514c53cc8d73d6b1078 100644
--- a/matlab/missing/mex/gensylv/gensylv.m
+++ b/matlab/missing/mex/gensylv/gensylv.m
@@ -1,5 +1,5 @@
-function [err, E] = gensylv(kron_prod,A,B,C0,D)
-%function [err, E] = gensylv(fake,A,B,C,D)
+function E = gensylv(kron_prod,A,B,C0,D)
+%function E = gensylv(fake,A,B,C,D)
 % Solves a Sylvester equation.
 %
 % INPUTS
@@ -19,7 +19,7 @@ function [err, E] = gensylv(kron_prod,A,B,C0,D)
 % SPECIAL REQUIREMENTS
 %   none.
 
-% Copyright (C) 1996-2017 Dynare Team
+% Copyright (C) 1996-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -42,4 +42,3 @@ end
 
 x0 = sylvester3(A,B,C,D);
 E  = sylvester3a(x0,A,B,C,D);
-err = 0;
\ No newline at end of file
diff --git a/matlab/missing/mex/kronecker/A_times_B_kronecker_C.m b/matlab/missing/mex/kronecker/A_times_B_kronecker_C.m
index 7a8d09498e5820f905ce68cbe509ca812e9dd9c4..db615145a1685a75771f830e3a2ed1b99c69209b 100644
--- a/matlab/missing/mex/kronecker/A_times_B_kronecker_C.m
+++ b/matlab/missing/mex/kronecker/A_times_B_kronecker_C.m
@@ -1,4 +1,4 @@
-function [D, err] = A_times_B_kronecker_C(A,B,C)
+function D = A_times_B_kronecker_C(A,B,C)
 
 %@info:
 %! @deftypefn {Function File} {[@var{D}, @var{err}] =} A_times_B_kronecker_C (@var{A},@var{B},@var{C})
@@ -41,7 +41,7 @@ function [D, err] = A_times_B_kronecker_C(A,B,C)
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 1996-2012 Dynare Team
+% Copyright (C) 1996-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -109,4 +109,4 @@ else
         D = A * kron(B,B);
     end
 end
-err = 0;
+
diff --git a/matlab/missing/mex/kronecker/sparse_hessian_times_B_kronecker_C.m b/matlab/missing/mex/kronecker/sparse_hessian_times_B_kronecker_C.m
index fd6059f840d9f0c2c9556b5abf06965f48437945..7857e2571e9167cb2634d1813608d87a93c6b777 100644
--- a/matlab/missing/mex/kronecker/sparse_hessian_times_B_kronecker_C.m
+++ b/matlab/missing/mex/kronecker/sparse_hessian_times_B_kronecker_C.m
@@ -1,7 +1,7 @@
-function [D, err] = sparse_hessian_times_B_kronecker_C(varargin)
+function D = sparse_hessian_times_B_kronecker_C(varargin)
 
 %@info:
-%! @deftypefn {Function File} {[@var{D}, @var{err}] =} sparse_hessian_times_B_kronecker_C (@var{A},@var{B},@var{C},@var{fake})
+%! @deftypefn {Function File} {@var{D} =} sparse_hessian_times_B_kronecker_C (@var{A},@var{B},@var{C},@var{fake})
 %! @anchor{kronecker/sparse_hessian_times_B_kronecker_C}
 %! @sp 1
 %! Computes A*kron(B,C) where A is hessian matrix in sparse format.
@@ -24,8 +24,6 @@ function [D, err] = sparse_hessian_times_B_kronecker_C(varargin)
 %! @table @ @var
 %! @item D
 %! mA*(nC*nB) or mA*(nB*nB) matrix of doubles.
-%! @item err
-%! Integer scalar equal to zero (if all goes well).
 %! @end table
 %! @sp 2
 %! @strong{Remarks}
@@ -45,7 +43,7 @@ function [D, err] = sparse_hessian_times_B_kronecker_C(varargin)
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 1996-2012 Dynare Team
+% Copyright (C) 1996-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -71,10 +69,9 @@ fake = varargin{nargin};
 
 switch nargin
   case 4
-    [D, fake] = A_times_B_kronecker_C(A,B,C,fake);
+    D = A_times_B_kronecker_C(A,B,C,fake);
   case 3
-    [D, fake] = A_times_B_kronecker_C(A,B,C);
+    D = A_times_B_kronecker_C(A,B,C);
   otherwise
     error('Two or Three input arguments required!')
 end
-err = 0;
\ No newline at end of file
diff --git a/matlab/missing/mex/mjdgges/mjdgges.m b/matlab/missing/mex/mjdgges/mjdgges.m
index 6bd345a570eaf0f4f37cd581be13a38b4bae164c..11d51d62f7709821b04a1a70fcd3fca381e6b262 100644
--- a/matlab/missing/mex/mjdgges/mjdgges.m
+++ b/matlab/missing/mex/mjdgges/mjdgges.m
@@ -1,4 +1,4 @@
-function [err, ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhreshold)
+function [ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhreshold)
 %
 % INPUTS
 %   e            [double] real square (n*n) matrix.
@@ -8,7 +8,6 @@ function [err, ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhr
 %                                  detecting eigenvalues too close to 0/0)
 %
 % OUTPUTS
-%   err          [integer] scalar: 1 indicates failure, 0 indicates success
 %   ss           [double] (n*n) quasi-triangular matrix.
 %   tt           [double] (n*n) quasi-triangular matrix.
 %   zz           [double] (n*n) orthogonal matrix.
@@ -19,7 +18,7 @@ function [err, ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhr
 % SPECIAL REQUIREMENTS
 %   none.
 
-% Copyright (C) 1996-2018 Dynare Team
+% Copyright (C) 1996-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -66,5 +65,3 @@ try
 catch
     info = 1; % Not as precise as lapack's info!
 end
-
-err = 0;
diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m
index 5696f772693ab9031acbd7d34c310cf85f69846f..70b39fe70772a8684e51c6076d17652261296278 100644
--- a/matlab/model_diagnostics.m
+++ b/matlab/model_diagnostics.m
@@ -16,7 +16,7 @@ function model_diagnostics(M,options,oo)
 %   none.
 %
 
-% Copyright (C) 1996-2018 Dynare Team
+% Copyright (C) 1996-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -107,11 +107,11 @@ exo = [oo.exo_steady_state; oo.exo_det_steady_state];
 for b=1:nb
     if options.bytecode
         if nb == 1
-            [chk, res, jacob] = bytecode(dr.ys, exo, M.params, dr.ys, 1, exo, ...
-                                         'evaluate', 'static');
+            [res, jacob] = bytecode(dr.ys, exo, M.params, dr.ys, 1, exo, ...
+                                    'evaluate', 'static');
         else
-            [chk, res, jacob] = bytecode(dr.ys, exo, M.params, dr.ys, 1, exo, ...
-                                         'evaluate', 'static',['block=' ...
+            [res, jacob] = bytecode(dr.ys, exo, M.params, dr.ys, 1, exo, ...
+                                    'evaluate', 'static',['block=' ...
                                 int2str(b)]);
         end
     else
@@ -204,8 +204,8 @@ z = repmat(dr.ys,1,klen);
 if ~options.block
     if options.order == 1
         if (options.bytecode)
-            [chck, ~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
-                                         M.params, dr.ys, 1);
+            [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
+                                   M.params, dr.ys, 1);
             jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd];
         else
             [~,jacobia_] = feval([M.fname '.dynamic'],z(iyr0),exo_simul, ...
@@ -213,8 +213,8 @@ if ~options.block
         end
     elseif options.order >= 2
         if (options.bytecode)
-            [chck, ~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
-                                         M.params, dr.ys, 1);
+            [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
+                                   M.params, dr.ys, 1);
             jacobia_ = [loc_dr.g1 loc_dr.g1_x];
         else
             [~,jacobia_,hessian1] = feval([M.fname '.dynamic'],z(iyr0),...
diff --git a/matlab/ms-sbvar/ms_compute_mdd.m b/matlab/ms-sbvar/ms_compute_mdd.m
index 499baf6a525b7495b6a3d6dd6e4116337d0d1dcc..43372a6e0a68851c1f6cade107cfb865513b9f3c 100644
--- a/matlab/ms-sbvar/ms_compute_mdd.m
+++ b/matlab/ms-sbvar/ms_compute_mdd.m
@@ -14,7 +14,7 @@ function [options_, oo_]=ms_compute_mdd(M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011-2013 Dynare Team
+% Copyright (C) 2011-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -51,8 +51,7 @@ if options_.ms.use_mean_center
 end
 
 % compute mdd
-[err] = ms_sbvar_command_line(opt);
-mexErrCheck('ms_compute_mdd',err);
+ms_sbvar_command_line(opt);
 
 % grab the muller/bridge log mdd from the output file
 mull_exp = 'Muller \w+\(\w+\) \= (\d+.\d+e\+\d+)';
diff --git a/matlab/ms-sbvar/ms_compute_probabilities.m b/matlab/ms-sbvar/ms_compute_probabilities.m
index 07b951f65ae5bf4209b318e89038bcd0aaae9b4b..9a5b91de891a82f4a46846f2e928b61fab03a381 100644
--- a/matlab/ms-sbvar/ms_compute_probabilities.m
+++ b/matlab/ms-sbvar/ms_compute_probabilities.m
@@ -14,7 +14,7 @@ function [options_, oo_]=ms_compute_probabilities(M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011 Dynare Team
+% Copyright (C) 2011-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -53,8 +53,7 @@ else
 end
 
 % compute probabilities
-[err] = ms_sbvar_command_line(opt);
-mexErrCheck('ms_compute_probabilities',err);
+ms_sbvar_command_line(opt);
 
 % now we want to plot the probabilities for each chain
 if ischar(prob_out_file)
diff --git a/matlab/ms-sbvar/ms_estimation.m b/matlab/ms-sbvar/ms_estimation.m
index 495a483eed3a00cbcc6e8821ddd10b4288353bbe..2a8510a2a440f95777fba32bb95bf95a658ab880 100644
--- a/matlab/ms-sbvar/ms_estimation.m
+++ b/matlab/ms-sbvar/ms_estimation.m
@@ -14,7 +14,7 @@ function [options_, oo_]=ms_estimation(M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011-2017 Dynare Team
+% Copyright (C) 2011-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -73,8 +73,7 @@ opt = [opt ' -random_tol_obj ' num2str(options_.ms.random_function_convergence_c
 opt = [opt ' -random_tol_parms ' num2str(options_.ms.random_parameter_convergence_criterion)];
 
 % estimation
-[err] = ms_sbvar_command_line(opt);
-mexErrCheck('ms_estimation', err);
+ms_sbvar_command_line(opt);
 
 [options_, oo_] = set_ms_estimation_file(options_.ms.output_file_tag, options_, oo_);
 end
diff --git a/matlab/ms-sbvar/ms_forecast.m b/matlab/ms-sbvar/ms_forecast.m
index e8949c9ae287d83700ad307f507752c3009094bb..b6e25aa847fc3f063b648422b7fff1f9463ce0e9 100644
--- a/matlab/ms-sbvar/ms_forecast.m
+++ b/matlab/ms-sbvar/ms_forecast.m
@@ -14,7 +14,7 @@ function [options_, oo_]=ms_forecast(M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011-2017 Dynare Team
+% Copyright (C) 2011-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -75,8 +75,7 @@ else
 end
 
 % forecast
-[err] = ms_sbvar_command_line(opt);
-mexErrCheck('ms_forecast',err);
+ms_sbvar_command_line(opt);
 
 % Plot Forecasts
 if options_.ms.regimes
diff --git a/matlab/ms-sbvar/ms_irf.m b/matlab/ms-sbvar/ms_irf.m
index 47d1a8490be13cde4d972f36c46cb66821f9aae0..d06a1149c0cead1e6abbea90fd6a8abd4a0ab1ab 100644
--- a/matlab/ms-sbvar/ms_irf.m
+++ b/matlab/ms-sbvar/ms_irf.m
@@ -15,7 +15,7 @@ function [options_, oo_]=ms_irf(varlist, M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011-2017 Dynare Team
+% Copyright (C) 2011-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -77,8 +77,7 @@ else
 end
 
 % irf
-[err] = ms_sbvar_command_line(opt);
-mexErrCheck('ms_irf',err);
+ms_sbvar_command_line(opt);
 
 % Plot IRFs
 if options_.ms.regimes
diff --git a/matlab/ms-sbvar/ms_sbvar_setup.m b/matlab/ms-sbvar/ms_sbvar_setup.m
index a195cc49137e0756fdde3c59bbde470808cd8ceb..062eaca728236b39f9b27b7a98c34734882d3f79 100644
--- a/matlab/ms-sbvar/ms_sbvar_setup.m
+++ b/matlab/ms-sbvar/ms_sbvar_setup.m
@@ -11,7 +11,7 @@ function ms_sbvar_setup(options_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2003-2017 Dynare Team
+% Copyright (C) 2003-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -428,6 +428,5 @@ fclose(fidForC);
 %======================================================================
 ms_write_markov_file(markov_file,options_)
 create_init_file = [matlab_filename,' ',markov_file,' ',options_.ms.file_tag];
-[err] = ms_sbvar_create_init_file(create_init_file);
-mexErrCheck('ms_sbvar_create_init_file',err);
+ms_sbvar_create_init_file(create_init_file);
 end
diff --git a/matlab/ms-sbvar/ms_simulation.m b/matlab/ms-sbvar/ms_simulation.m
index f7832e035f372de8a683b0be133ff7c74790cc65..d311733b80c9712197c95c09d1f39a81d6bca2a3 100644
--- a/matlab/ms-sbvar/ms_simulation.m
+++ b/matlab/ms-sbvar/ms_simulation.m
@@ -14,7 +14,7 @@ function [options_, oo_]=ms_simulation(M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011-2013 Dynare Team
+% Copyright (C) 2011-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -50,6 +50,5 @@ if options_.ms.save_draws
 end
 
 % simulation
-[err] = ms_sbvar_command_line(opt);
-mexErrCheck('ms_simulation',err);
+ms_sbvar_command_line(opt);
 end
diff --git a/matlab/ms-sbvar/ms_variance_decomposition.m b/matlab/ms-sbvar/ms_variance_decomposition.m
index 1fe73ddbc4edab37d68b145527974818635e158b..746facf2ca322051f28c3929161dff391ff85417 100644
--- a/matlab/ms-sbvar/ms_variance_decomposition.m
+++ b/matlab/ms-sbvar/ms_variance_decomposition.m
@@ -14,7 +14,7 @@ function [options_, oo_]=ms_variance_decomposition(M_, options_, oo_)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2011-2017 Dynare Team
+% Copyright (C) 2011-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -78,8 +78,7 @@ if options_.ms.error_bands
 end
 
 % variance_decomposition
-[err] = ms_sbvar_command_line(opt);
-mexErrCheck('ms_variance_decomposition',err);
+ms_sbvar_command_line(opt);
 
 if options_.ms.regime || options_.ms.regimes
     outfile = [outfile 'regime_'];
diff --git a/matlab/partial_information/add_auxiliary_variables_to_steadystate.m b/matlab/partial_information/add_auxiliary_variables_to_steadystate.m
index e1ffe4f196e9eefc9c0bf3d40a1d2d6f2ef23094..35bd84bb33e794c93508f4d5d92074dd0596a52d 100644
--- a/matlab/partial_information/add_auxiliary_variables_to_steadystate.m
+++ b/matlab/partial_information/add_auxiliary_variables_to_steadystate.m
@@ -2,7 +2,7 @@ function ys1 = add_auxiliary_variables_to_steadystate(ys,aux_vars,fname, ...
                                                   exo_steady_state, exo_det_steady_state,params, byte_code)
 % Add auxiliary variables to the steady state vector
 
-% Copyright (C) 2009-2017 Dynare Team
+% Copyright (C) 2009-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -23,9 +23,9 @@ ys1 = [ys;zeros(n,1)];
 
 for i=1:n+1
     if byte_code
-        [info, res] = bytecode('static','evaluate',ys1,...
-                               [exo_steady_state; ...
-                            exo_det_steady_state],params);
+        res = bytecode('static','evaluate',ys1,...
+                       [exo_steady_state; ...
+                        exo_det_steady_state],params);
     else
         res = feval([fname '.static'],ys1,...
                     [exo_steady_state; ...
diff --git a/matlab/perfect-foresight-models/det_cond_forecast.m b/matlab/perfect-foresight-models/det_cond_forecast.m
index 98875fb91757d54442cab2f511e9be3f5d97ca62..cdcbe24d77132486930870151c6efedc49c7bfb4 100644
--- a/matlab/perfect-foresight-models/det_cond_forecast.m
+++ b/matlab/perfect-foresight-models/det_cond_forecast.m
@@ -12,7 +12,7 @@ function data_set = det_cond_forecast(varargin)
 %  dataset                [dseries]     Returns a dseries containing the forecasted endgenous variables and shocks
 %
 
-% Copyright (C) 2013-2018 Dynare Team
+% Copyright (C) 2013-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -158,13 +158,12 @@ else
             if options_.bytecode
                 save_options_dynatol_f = options_.dynatol.f;
                 options_.dynatol.f = 1e-7;
-                [Info, endo, exo] = bytecode('extended_path', plan, oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, options_.periods);
+                [endo, exo] = bytecode('extended_path', plan, oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, options_.periods);
                 options_.dynatol.f = save_options_dynatol_f;
 
-                if Info == 0
-                    oo_.endo_simul = endo;
-                    oo_.exo_simul = exo;
-                end
+                oo_.endo_simul = endo;
+                oo_.exo_simul = exo;
+
                 endo = endo';
                 endo_l = size(endo(1+M_.maximum_lag:end,:),1);
                 jrng = dates(plan.date(1)):dates(plan.date(1)+endo_l);
@@ -473,14 +472,12 @@ if pf && ~surprise
                 end
                 data1 = M_;
                 if (options_.bytecode)
-                    [chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
+                    [zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
                 else
                     [zz, g1b] = feval([M_.fname '.dynamic'], z', zx, M_.params, oo_.steady_state, k);
                     data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end);
                     data1.g1 = g1b(:,1 : end - M_.exo_nbr);
-                    chck = 0;
                 end
-                mexErrCheck('bytecode', chck);
             end
             if k == 1
                 g1(1:M_.endo_nbr,-M_.endo_nbr + [cur_indx lead_indx]) = data1.g1(:,M_.nspred + 1:end);
@@ -745,14 +742,12 @@ else
                     end
                     data1 = M_;
                     if (options_.bytecode)
-                        [chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
+                        [zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
                     else
                         [zz, g1b] = feval([M_.fname '.dynamic'], z', zx, M_.params, oo_.steady_state, k);
                         data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end);
                         data1.g1 = g1b(:,1 : end - M_.exo_nbr);
-                        chck = 0;
                     end
-                    mexErrCheck('bytecode', chck);
                 end
                 if k == 1
                     g1(1:M_.endo_nbr,-M_.endo_nbr + [cur_indx lead_indx]) = data1.g1(:,M_.nspred + 1:end);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 85b420cd4bca22b088d88b9d833c511b3285ccd3..267545d69247cef70c66adaee222edf8068040ad 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -11,7 +11,7 @@ function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_)
 % - oo_                 [struct] contains results
 % - maxerror            [double] contains the maximum absolute error
 
-% Copyright (C) 2015-2019 Dynare Team
+% Copyright (C) 2015-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -50,18 +50,14 @@ end
 if options_.block
     if options_.bytecode
         try
-            [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
-        catch
-            info = 1;
-        end
-        if info
-            oo_.deterministic_simulation.status = false;
-        else
-            oo_.endo_simul = tmp;
+            oo_.endo_simul = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
             oo_.deterministic_simulation.status = true;
-        end
-        if options_.no_homotopy
-            mexErrCheck('bytecode', info);
+        catch ME
+            disp(ME.message)
+            if options_.no_homotopy
+                error('Error in bytecode')
+            end
+            oo_.deterministic_simulation.status = false;
         end
     else
         oo_ = feval([M_.fname '.dynamic'], options_, M_, oo_);
@@ -69,18 +65,14 @@ if options_.block
 else
     if options_.bytecode
         try
-            [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods);
-        catch
-            info = 1;
-        end
-        if info
-            oo_.deterministic_simulation.status = false;
-        else
-            oo_.endo_simul = tmp;
+            oo_.endo_simul = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods);
             oo_.deterministic_simulation.status = true;
-        end
-        if options_.no_homotopy
-            mexErrCheck('bytecode', info);
+        catch ME
+            disp(ME.message)
+            if options_.no_homotopy
+                error('Error in bytecode')
+            end
+            oo_.deterministic_simulation.status = false;
         end
     else
         if M_.maximum_endo_lead == 0 % Purely backward model
@@ -128,7 +120,7 @@ if nargout>1
         maxerror = oo_.deterministic_simulation.error;
     else
         if options_.bytecode
-            [~, residuals]= bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
+            residuals = bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
         else
             if M_.maximum_lag > 0
                 y0 = oo_.endo_simul(:, M_.maximum_lag);
diff --git a/matlab/resid.m b/matlab/resid.m
index 716c9c5b36f721d88ed92719ed3647d788a9d10f..1f3743a072e7c590931eb1e9f0986414b672738a 100644
--- a/matlab/resid.m
+++ b/matlab/resid.m
@@ -12,7 +12,7 @@ function z = resid(junk)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2001-2018 Dynare Team
+% Copyright (C) 2001-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -79,8 +79,7 @@ if options_.block && ~options_.bytecode
         z(idx) = r;
     end
 elseif options_.bytecode
-    [check, z] = bytecode('evaluate','static');
-    mexErrCheck('bytecode', check);
+    z = bytecode('evaluate','static');
 else
     z = feval([M_.fname '.static'],...
               oo_.steady_state,...
diff --git a/matlab/simult_.m b/matlab/simult_.m
index e3ec1ec31724d2e9d4bd9816a2f6b38c2a1e9af2..508af26e435fe890db154dd1a1f9733c29b46911 100644
--- a/matlab/simult_.m
+++ b/matlab/simult_.m
@@ -17,7 +17,7 @@ function y_=simult_(M_,options_,y0,dr,ex_,iorder)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2001-2019 Dynare Team
+% Copyright (C) 2001-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -53,10 +53,9 @@ end
 
 if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
     ex_ = [zeros(M_.maximum_lag,M_.exo_nbr); ex_];
-    [err, y_] = dynare_simul_(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,exo_nbr, ...
-                              y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed, ...
-                              dr.ys(dr.order_var),dr);
-    mexErrCheck('dynare_simul_', err);
+    y_ = dynare_simul_(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,exo_nbr, ...
+                       y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed, ...
+                       dr.ys(dr.order_var),dr);
     y_(dr.order_var,:) = y_;
 else
     if options_.block
@@ -98,12 +97,9 @@ else
                 yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
                 yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
                 epsilon = ex_(i-1,:)';
-                [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1);
-                mexErrCheck('A_times_B_kronecker_C', err);
-                [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
-                mexErrCheck('A_times_B_kronecker_C', err);
-                [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon);
-                mexErrCheck('A_times_B_kronecker_C', err);
+                abcOut1 = A_times_B_kronecker_C(.5*dr.ghxx,yhat1);
+                abcOut2 = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
+                abcOut3 = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon);
                 y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
                     + abcOut1 + abcOut2 + abcOut3;
                 y__(order_var) = dr.ys(order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
@@ -112,12 +108,9 @@ else
             for i = 2:iter+M_.maximum_lag
                 yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
                 epsilon = ex_(i-1,:)';
-                [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat);
-                mexErrCheck('A_times_B_kronecker_C', err);
-                [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
-                mexErrCheck('A_times_B_kronecker_C', err);
-                [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon);
-                mexErrCheck('A_times_B_kronecker_C', err);
+                abcOut1 = A_times_B_kronecker_C(.5*dr.ghxx,yhat);
+                abcOut2 = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
+                abcOut3 = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon);
                 y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
                     + abcOut1 + abcOut2 + abcOut3;
             end
@@ -150,30 +143,21 @@ else
             u = ex_(i-1,:)';
             %construct terms of order 2 from second order part, based
             %on linear component yhat1
-            [gyy, err] = A_times_B_kronecker_C(ghxx,yhat1);
-            mexErrCheck('A_times_B_kronecker_C', err);
-            [guu, err] = A_times_B_kronecker_C(ghuu,u);
-            mexErrCheck('A_times_B_kronecker_C', err);
-            [gyu, err] = A_times_B_kronecker_C(ghxu,yhat1,u);
-            mexErrCheck('A_times_B_kronecker_C', err);
+            gyy = A_times_B_kronecker_C(ghxx,yhat1);
+            guu = A_times_B_kronecker_C(ghuu,u);
+            gyu = A_times_B_kronecker_C(ghxu,yhat1,u);
             %construct terms of order 3 from second order part, based
             %on order 2 component yhat2
-            [gyy12, err] = A_times_B_kronecker_C(ghxx,yhat1,yhat2);
-            mexErrCheck('A_times_B_kronecker_C', err);
-            [gy2u, err] = A_times_B_kronecker_C(ghxu,yhat2,u);
-            mexErrCheck('A_times_B_kronecker_C', err);
+            gyy12 = A_times_B_kronecker_C(ghxx,yhat1,yhat2);
+            gy2u = A_times_B_kronecker_C(ghxu,yhat2,u);
             %construct terms of order 3, all based on first order component yhat1
             y2a = kron(yhat1,yhat1);
-            [gyyy, err] = A_times_B_kronecker_C(ghxxx,y2a,yhat1);
-            mexErrCheck('A_times_B_kronecker_C', err);
+            gyyy = A_times_B_kronecker_C(ghxxx,y2a,yhat1);
             u2a = kron(u,u);
-            [guuu, err] = A_times_B_kronecker_C(ghuuu,u2a,u);
-            mexErrCheck('A_times_B_kronecker_C', err);
+            guuu = A_times_B_kronecker_C(ghuuu,u2a,u);
             yu = kron(yhat1,u);
-            [gyyu, err] = A_times_B_kronecker_C(ghxxu,yhat1,yu);
-            mexErrCheck('A_times_B_kronecker_C', err);
-            [gyuu, err] = A_times_B_kronecker_C(ghxuu,yu,u);
-            mexErrCheck('A_times_B_kronecker_C', err);
+            gyyu = A_times_B_kronecker_C(ghxxu,yhat1,yu);
+            gyuu = A_times_B_kronecker_C(ghxuu,yu,u);
             %add all terms of order 3, linear component based on third
             %order yhat3
             yhat3 = ghx*yhat3 +gyy12 ... % prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
diff --git a/matlab/solve_perfect_foresight_model.m b/matlab/solve_perfect_foresight_model.m
index cedf3e0604f0419d58e95514ebd1085c5a2e93b1..93c5fb26fbcaae6e6d9614e9c82ea5151cf07e2f 100644
--- a/matlab/solve_perfect_foresight_model.m
+++ b/matlab/solve_perfect_foresight_model.m
@@ -1,6 +1,6 @@
 function [flag,endo_simul,err] = solve_perfect_foresight_model(endo_simul,exo_simul,pfm)
 
-% Copyright (C) 2012-2017 Dynare Team
+% Copyright (C) 2012-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,7 +33,13 @@ if pfm.verbose
 end
 
 if pfm.use_bytecode
-    [flag, endo_simul]=bytecode(Y, exo_simul, pfm.params);
+    try
+        endo_simul=bytecode(Y, exo_simul, pfm.params);
+        flag = 0;
+    catch ME
+        disp(ME.message);
+        flag = 1;
+    end
     return
 end
 
@@ -120,4 +126,4 @@ if nan_flag
 end
 if pfm.verbose
     disp (['-----------------------------------------------------']) ;
-end
\ No newline at end of file
+end
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index 5aa9920b410c7032925704597380eabbeff6a2bb..87d88bacf128e8f6cbebe51320344cced1a7d693 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -24,7 +24,7 @@ function [dr, info] = stochastic_solvers(dr, task, M_, options_, oo_)
 %                                 info=6 -> The jacobian matrix evaluated at the steady state is complex.
 %                                 info=9 -> k_order_pert was unable to compute the solution
 
-% Copyright (C) 1996-2019 Dynare Team
+% Copyright (C) 1996-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -111,8 +111,8 @@ it_ = M_.maximum_lag + 1;
 z = repmat(dr.ys,1,klen);
 if local_order == 1
     if (options_.bytecode)
-        [chck, ~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
-                                     M_.params, dr.ys, 1);
+        [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
+                               M_.params, dr.ys, 1);
         jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd];
     else
         [~,jacobia_] = feval([M_.fname '.dynamic'],z(iyr0),exo_simul, ...
@@ -120,8 +120,8 @@ if local_order == 1
     end
 elseif local_order == 2
     if (options_.bytecode)
-        [chck, ~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
-                                     M_.params, dr.ys, 1);
+        [~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
+                               M_.params, dr.ys, 1);
         jacobia_ = [loc_dr.g1 loc_dr.g1_x];
     else
         [~,jacobia_,hessian1] = feval([M_.fname '.dynamic'],z(iyr0),...
diff --git a/mex/sources/block_kalman_filter/block_kalman_filter.cc b/mex/sources/block_kalman_filter/block_kalman_filter.cc
index 0f484b04ffc2c60bf87a40ea828d28e9d14e6114..0ad1f60c80f17b9549d97dc4593fcc85ee644d4a 100644
--- a/mex/sources/block_kalman_filter/block_kalman_filter.cc
+++ b/mex/sources/block_kalman_filter/block_kalman_filter.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2019 Dynare Team
+ * Copyright © 2007-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -143,10 +143,10 @@ det(const double *F, int dim, const lapack_int *ipiv)
 
 BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
-  if (nlhs > 3)
-    DYN_MEX_FUNC_ERR_MSG_TXT("block_kalman_filter provides at most 3 output argument.");
+  if (nlhs > 2)
+    mexErrMsgTxt("block_kalman_filter provides at most 2 output argument.");
   if (nrhs != 13 && nrhs != 16)
-    DYN_MEX_FUNC_ERR_MSG_TXT("block_kalman_filter requires exactly \n  13 input arguments for standard Kalman filter \nor\n  16 input arguments for missing observations Kalman filter.");
+    mexErrMsgTxt("block_kalman_filter requires exactly \n  13 input arguments for standard Kalman filter \nor\n  16 input arguments for missing observations Kalman filter.");
   if (nrhs == 16)
     missing_observations = true;
   else
@@ -154,13 +154,13 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
   if (missing_observations)
     {
       if (!mxIsCell(prhs[0]))
-        DYN_MEX_FUNC_ERR_MSG_TXT("the first input argument of block_missing_observations_kalman_filter must be a Cell Array.");
+        mexErrMsgTxt("the first input argument of block_missing_observations_kalman_filter must be a Cell Array.");
       pdata_index = prhs[0];
       if (!mxIsDouble(prhs[1]))
-        DYN_MEX_FUNC_ERR_MSG_TXT("the second input argument of block_missing_observations_kalman_filter must be a scalar.");
+        mexErrMsgTxt("the second input argument of block_missing_observations_kalman_filter must be a scalar.");
       number_of_observations = ceil(mxGetScalar(prhs[1]));
       if (!mxIsDouble(prhs[2]))
-        DYN_MEX_FUNC_ERR_MSG_TXT("the third input argument of block_missing_observations_kalman_filter must be a scalar.");
+        mexErrMsgTxt("the third input argument of block_missing_observations_kalman_filter must be a scalar.");
       no_more_missing_observations = ceil(mxGetScalar(prhs[2]));
       pT = mxDuplicateArray(prhs[3]);
       pR = mxDuplicateArray(prhs[4]);
@@ -218,7 +218,7 @@ BlockKalmanFilter::BlockKalmanFilter(int nlhs, mxArray *plhs[], int nrhs, const
 
   if (missing_observations)
     if (mxGetNumberOfElements(pdata_index) != static_cast<unsigned int>(smpl))
-      DYN_MEX_FUNC_ERR_MSG_TXT("the number of element in the cell array passed to block_missing_observation_kalman_filter as first argument has to be equal to the smpl size");
+      mexErrMsgTxt("the number of element in the cell array passed to block_missing_observation_kalman_filter as first argument has to be equal to the smpl size");
 
   i_nz_state_var = std::make_unique<int[]>(n);
   for (int i = 0; i < n; i++)
@@ -436,7 +436,7 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[])
           {
             mexPrintf("error: F singular\n");
             LIK = Inf;
-            if (nlhs == 3)
+            if (nlhs == 2)
               for (int i = t; i < smpl; i++)
                 lik[i] = Inf;
             // info = 0
@@ -817,17 +817,15 @@ BlockKalmanFilter::block_kalman_filter(int nlhs, mxArray *plhs[])
 void
 BlockKalmanFilter::return_results_and_clean(int nlhs, mxArray *plhs[])
 {
-  plhs[0] = mxCreateDoubleScalar(0);
-
-  if (nlhs >= 2)
+  if (nlhs >= 1)
     {
-      plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
-      double *pind = mxGetPr(plhs[1]);
+      plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+      double *pind = mxGetPr(plhs[0]);
       pind[0] = LIK;
     }
 
-  if (nlhs == 3)
-    plhs[2] = plik;
+  if (nlhs == 2)
+    plhs[1] = plik;
   else
     mxDestroyArray(plik);
 
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index ee81f0843cd5982f57fd12b36656fe51eed12b70..31efcd311c6f683b0cad3973f98203ba2f62abb4 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2017 Dynare Team
+ * Copyright © 2007-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -21,37 +21,6 @@
 #include "ErrorHandling.hh"
 #include <ctime>
 #include <math.h>
-#ifdef DYN_MEX_FUNC_ERR_MSG_TXT
-# undef DYN_MEX_FUNC_ERR_MSG_TXT
-#endif // DYN_MEX_FUNC_ERR_MSG_TXT
-
-#define DYN_MEX_FUNC_ERR_MSG_TXT(str)                                   \
-  do {                                                                  \
-    mexPrintf("%s\n", str);                                             \
-    if (nlhs > 0)                                                       \
-      {                                                                 \
-        plhs[0] = mxCreateDoubleScalar(1);                              \
-        if (nlhs > 1)                                                   \
-          {                                                             \
-            double *pind;                                               \
-            plhs[1] = mxCreateDoubleMatrix(int (row_y), int (col_y), mxREAL); \
-            pind = mxGetPr(plhs[1]);                                    \
-            if (evaluate)                                               \
-              {                                                         \
-                for (unsigned int i = 0; i < row_y*col_y; i++)          \
-                  pind[i] = 0;                                          \
-              }                                                         \
-            else                                                        \
-              {                                                         \
-                for (unsigned int i = 0; i < row_y*col_y; i++)          \
-                  pind[i] = yd[i];                                      \
-              }                                                         \
-            for (int i = 2; i < nlhs; i++)                              \
-              plhs[i] = mxCreateDoubleScalar(1);                        \
-          }                                                             \
-      }                                                                 \
-    return;                                                             \
-  } while (0)
 
 #ifdef DEBUG_EX
 
@@ -492,13 +461,13 @@ main(int nrhs, const char *prhs[])
     }
   catch (GeneralExceptionHandling &feh)
     {
-      DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
+      mexErrMsgTxt(feh.GetErrorMsg().c_str());
     }
   if (!count_array_argument)
     {
       int field = mxGetFieldNumber(M_, "params");
       if (field < 0)
-        DYN_MEX_FUNC_ERR_MSG_TXT("params is not a field of M_");
+        mexErrMsgTxt("params is not a field of M_");
       params = mxGetPr(mxGetFieldByNumber(M_, 0, field));
     }
 
@@ -510,14 +479,14 @@ main(int nrhs, const char *prhs[])
       if (extended_path_struct == NULL)
         {
           string tmp = "The 'extended_path' option must be followed by the extended_path descriptor";
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       mxArray *date_str = mxGetField(extended_path_struct, 0, "date_str");
       if (date_str == NULL)
         {
           string tmp = "date_str";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       int nb_periods = mxGetM(date_str) * mxGetN(date_str);
 
@@ -526,28 +495,28 @@ main(int nrhs, const char *prhs[])
         {
           string tmp = "constrained_vars_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       mxArray *constrained_paths_ = mxGetField(extended_path_struct, 0, "constrained_paths_");
       if (constrained_paths_ == NULL)
         {
           string tmp = "constrained_paths_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       mxArray *constrained_int_date_ = mxGetField(extended_path_struct, 0, "constrained_int_date_");
       if (constrained_int_date_ == NULL)
         {
           string tmp = "constrained_int_date_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       mxArray *constrained_perfect_foresight_ = mxGetField(extended_path_struct, 0, "constrained_perfect_foresight_");
       if (constrained_perfect_foresight_ == NULL)
         {
           string tmp = "constrained_perfect_foresight_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
 
       mxArray *shock_var_ = mxGetField(extended_path_struct, 0, "shock_vars_");
@@ -555,28 +524,28 @@ main(int nrhs, const char *prhs[])
         {
           string tmp = "shock_vars_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       mxArray *shock_paths_ = mxGetField(extended_path_struct, 0, "shock_paths_");
       if (shock_paths_ == NULL)
         {
           string tmp = "shock_paths_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       mxArray *shock_int_date_ = mxGetField(extended_path_struct, 0, "shock_int_date_");
       if (shock_int_date_ == NULL)
         {
           string tmp = "shock_int_date_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       mxArray *shock_str_date_ = mxGetField(extended_path_struct, 0, "shock_str_date_");
       if (shock_str_date_ == NULL)
         {
           string tmp = "shock_str_date_";
           tmp.insert(0, "The extended_path description structure does not contain the member: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       int nb_constrained = mxGetM(constrained_vars_) * mxGetN(constrained_vars_);
       int nb_controlled = 0;
@@ -588,7 +557,7 @@ main(int nrhs, const char *prhs[])
           nb_controlled = mxGetM(controlled_varexo) * mxGetN(controlled_varexo);
           if (nb_controlled != nb_constrained)
             {
-              DYN_MEX_FUNC_ERR_MSG_TXT("The number of exogenized variables and the number of exogenous controlled variables should be equal.");
+              mexErrMsgTxt("The number of exogenized variables and the number of exogenous controlled variables should be equal.");
             }
         }
       double *controlled_varexo_value = NULL;
@@ -637,7 +606,7 @@ main(int nrhs, const char *prhs[])
               string tmp1 = oss.str();
               tmp.append(tmp1);
               tmp.append(")");
-              DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+              mexErrMsgTxt(tmp.c_str());
             }
           (sconditional_extended_path[i]).per_value.resize(nb_local_periods);
           (sconditional_extended_path[i]).value.resize(nb_periods);
@@ -680,7 +649,7 @@ main(int nrhs, const char *prhs[])
               string tmp1 = oss.str();
               tmp.append(tmp1);
               tmp.append(")");
-              DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+              mexErrMsgTxt(tmp.c_str());
             }
           (sextended_path[i]).per_value.resize(nb_local_periods);
           (sextended_path[i]).value.resize(nb_periods);
@@ -702,7 +671,7 @@ main(int nrhs, const char *prhs[])
           if (info)
             {
               string tmp = "Can not allocated memory to store the date_str in the extended path descriptor";
-              DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+              mexErrMsgTxt(tmp.c_str());
             }
           dates.push_back(string(buf)); //string(Dates[i]);
           mxFree(buf);
@@ -715,7 +684,7 @@ main(int nrhs, const char *prhs[])
         {
           string tmp = plan;
           tmp.insert(0, "Can't find the plan: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       size_t n_plan = mxGetN(plan_struct);
       splan.resize(n_plan);
@@ -738,7 +707,7 @@ main(int nrhs, const char *prhs[])
                   string tmp = name;
                   tmp.insert(0, "the variable '");
                   tmp.append("'  defined as var in plan is not an exogenous or a deterministic exogenous\n");
-                  DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+                  mexErrMsgTxt(tmp.c_str());
                 }
             }
           tmp = mxGetField(plan_struct, i, "var");
@@ -756,7 +725,7 @@ main(int nrhs, const char *prhs[])
                   string tmp = name;
                   tmp.insert(0, "the variable '");
                   tmp.append("'  defined as exo in plan is not an endogenous variable\n");
-                  DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+                  mexErrMsgTxt(tmp.c_str());
                 }
             }
           tmp = mxGetField(plan_struct, i, "per_value");
@@ -793,7 +762,7 @@ main(int nrhs, const char *prhs[])
         {
           string tmp = pfplan;
           tmp.insert(0, "Can't find the pfplan: ");
-          DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+          mexErrMsgTxt(tmp.c_str());
         }
       size_t n_plan = mxGetN(pfplan_struct);
       spfplan.resize(n_plan);
@@ -816,7 +785,7 @@ main(int nrhs, const char *prhs[])
                   string tmp = name;
                   tmp.insert(0, "the variable '");
                   tmp.append("' defined as var in pfplan is not an exogenous or a deterministic exogenous\n");
-                  DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+                  mexErrMsgTxt(tmp.c_str());
                 }
             }
           tmp = mxGetField(pfplan_struct, i, "exo");
@@ -834,7 +803,7 @@ main(int nrhs, const char *prhs[])
                   string tmp = name;
                   tmp.insert(0, "the variable '");
                   tmp.append("' defined as exo in pfplan  is not an endogenous variable\n");
-                  DYN_MEX_FUNC_ERR_MSG_TXT(tmp.c_str());
+                  mexErrMsgTxt(tmp.c_str());
                 }
             }
           tmp = mxGetField(pfplan_struct, i, "per_value");
@@ -866,20 +835,20 @@ main(int nrhs, const char *prhs[])
 
   int field_steady_state = mxGetFieldNumber(oo_, "steady_state");
   if (field_steady_state < 0)
-    DYN_MEX_FUNC_ERR_MSG_TXT("steady_state is not a field of oo_");
+    mexErrMsgTxt("steady_state is not a field of oo_");
   int field_exo_steady_state = mxGetFieldNumber(oo_, "exo_steady_state");
   if (field_exo_steady_state < 0)
-    DYN_MEX_FUNC_ERR_MSG_TXT("exo_steady_state is not a field of oo_");
+    mexErrMsgTxt("exo_steady_state is not a field of oo_");
 
   if (!steady_state)
     {
       int field_endo_simul = mxGetFieldNumber(oo_, "endo_simul");
       if (field_endo_simul < 0)
-        DYN_MEX_FUNC_ERR_MSG_TXT("endo_simul is not a field of oo_");
+        mexErrMsgTxt("endo_simul is not a field of oo_");
 
       int field_exo_simul = mxGetFieldNumber(oo_, "exo_simul");
       if (field_exo_simul < 0)
-        DYN_MEX_FUNC_ERR_MSG_TXT("exo_simul is not a field of oo_");
+        mexErrMsgTxt("exo_simul is not a field of oo_");
 
       if (!count_array_argument)
         {
@@ -897,17 +866,17 @@ main(int nrhs, const char *prhs[])
       if (field >= 0)
         y_kmin = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("maximum_lag is not a field of M_");
+        mexErrMsgTxt("maximum_lag is not a field of M_");
       field = mxGetFieldNumber(M_, "maximum_lead");
       if (field >= 0)
         y_kmax = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field)))));
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("maximum_lead is not a field of M_");
+        mexErrMsgTxt("maximum_lead is not a field of M_");
       field = mxGetFieldNumber(M_, "maximum_endo_lag");
       if (field >= 0)
         y_decal = max(0, y_kmin-int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, field))))));
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("maximum_endo_lag is not a field of M_");
+        mexErrMsgTxt("maximum_endo_lag is not a field of M_");
 
       if (!count_array_argument)
         {
@@ -915,7 +884,7 @@ main(int nrhs, const char *prhs[])
           if (field >= 0)
             periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, field)))));
           else
-            DYN_MEX_FUNC_ERR_MSG_TXT("options_ is not a field of options_");
+            mexErrMsgTxt("options_ is not a field of options_");
         }
 
       if (!steady_yd)
@@ -948,7 +917,7 @@ main(int nrhs, const char *prhs[])
   if (field >= 0)
     verbose = int (*mxGetPr((mxGetFieldByNumber(options_, 0, field))));
   else
-    DYN_MEX_FUNC_ERR_MSG_TXT("verbosity is not a field of options_");
+    mexErrMsgTxt("verbosity is not a field of options_");
   if (verbose)
     print_it = true;
   if (!steady_state)
@@ -961,34 +930,34 @@ main(int nrhs, const char *prhs[])
   else
     {
       if (!steady_state)
-        DYN_MEX_FUNC_ERR_MSG_TXT("simul is not a field of options_");
+        mexErrMsgTxt("simul is not a field of options_");
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("steady is not a field of options_");
+        mexErrMsgTxt("steady is not a field of options_");
     }
   field = mxGetFieldNumber(temporaryfield, "maxit");
   if (field < 0)
     {
       if (!steady_state)
-        DYN_MEX_FUNC_ERR_MSG_TXT("maxit is not a field of options_.simul");
+        mexErrMsgTxt("maxit is not a field of options_.simul");
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("maxit is not a field of options_.steady");
+        mexErrMsgTxt("maxit is not a field of options_.steady");
     }
   int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(temporaryfield, 0, field)))));
   field = mxGetFieldNumber(options_, "slowc");
   if (field < 0)
-    DYN_MEX_FUNC_ERR_MSG_TXT("slows is not a field of options_");
+    mexErrMsgTxt("slows is not a field of options_");
   double slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
   field = mxGetFieldNumber(options_, "markowitz");
   if (field < 0)
-    DYN_MEX_FUNC_ERR_MSG_TXT("markowitz is not a field of options_");
+    mexErrMsgTxt("markowitz is not a field of options_");
   double markowitz_c = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
   field = mxGetFieldNumber(options_, "minimal_solving_periods");
   if (field < 0)
-    DYN_MEX_FUNC_ERR_MSG_TXT("minimal_solving_periods is not a field of options_");
+    mexErrMsgTxt("minimal_solving_periods is not a field of options_");
   int minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
   field = mxGetFieldNumber(options_, "stack_solve_algo");
   if (field < 0)
-    DYN_MEX_FUNC_ERR_MSG_TXT("stack_solve_algo is not a field of options_");
+    mexErrMsgTxt("stack_solve_algo is not a field of options_");
   int stack_solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
   int solve_algo;
   double solve_tolf;
@@ -999,12 +968,12 @@ main(int nrhs, const char *prhs[])
       if (field >= 0)
         solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, field))));
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("solve_algo is not a field of options_");
+        mexErrMsgTxt("solve_algo is not a field of options_");
       field = mxGetFieldNumber(options_, "solve_tolf");
       if (field >= 0)
         solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, field)));
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("solve_tolf is not a field of options_");
+        mexErrMsgTxt("solve_tolf is not a field of options_");
     }
   else
     {
@@ -1014,19 +983,19 @@ main(int nrhs, const char *prhs[])
       if (field >= 0)
         dynatol = mxGetFieldByNumber(options_, 0, field);
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("dynatol is not a field of options_");
+        mexErrMsgTxt("dynatol is not a field of options_");
       field = mxGetFieldNumber(dynatol, "f");
       if (field >= 0)
         solve_tolf = *mxGetPr((mxGetFieldByNumber(dynatol, 0, field)));
       else
-        DYN_MEX_FUNC_ERR_MSG_TXT("f is not a field of options_.dynatol");
+        mexErrMsgTxt("f is not a field of options_.dynatol");
     }
   field = mxGetFieldNumber(M_, "fname");
   mxArray *mxa;
   if (field >= 0)
     mxa = mxGetFieldByNumber(M_, 0, field);
   else
-    DYN_MEX_FUNC_ERR_MSG_TXT("fname is not a field of M_");
+    mexErrMsgTxt("fname is not a field of M_");
   size_t buflen = mxGetM(mxa) * mxGetN(mxa) + 1;
   char *fname;
   fname = static_cast<char *>(mxCalloc(buflen+1, sizeof(char)));
@@ -1044,11 +1013,11 @@ main(int nrhs, const char *prhs[])
     }
   catch (GeneralExceptionHandling &feh)
     {
-      DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
+      mexErrMsgTxt(feh.GetErrorMsg().c_str());
     }
 #else
   if (stack_solve_algo == 7 && !steady_state)
-    DYN_MEX_FUNC_ERR_MSG_TXT("bytecode has not been compiled with CUDA option. Bytecode Can't use options_.stack_solve_algo=7\n");
+    mexErrMsgTxt("bytecode has not been compiled with CUDA option. Bytecode Can't use options_.stack_solve_algo=7\n");
 #endif
   size_t size_of_direction = col_y*row_y*sizeof(double);
   double *y = static_cast<double *>(mxMalloc(size_of_direction));
@@ -1085,7 +1054,6 @@ main(int nrhs, const char *prhs[])
   mxFree(fname);
   int nb_blocks = 0;
   double *pind;
-  bool no_error = true;
 
   if (extended_path)
     {
@@ -1095,7 +1063,7 @@ main(int nrhs, const char *prhs[])
         }
       catch (GeneralExceptionHandling &feh)
         {
-          DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
+          mexErrMsgTxt(feh.GetErrorMsg().c_str());
         }
     }
   else
@@ -1106,7 +1074,7 @@ main(int nrhs, const char *prhs[])
         }
       catch (GeneralExceptionHandling &feh)
         {
-          DYN_MEX_FUNC_ERR_MSG_TXT(feh.GetErrorMsg().c_str());
+          mexErrMsgTxt(feh.GetErrorMsg().c_str());
         }
     }
 
@@ -1116,42 +1084,21 @@ main(int nrhs, const char *prhs[])
 #endif
 
   clock_t t1 = clock();
-  if (!steady_state && !evaluate && no_error && print)
+  if (!steady_state && !evaluate && print)
     mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
 #ifndef DEBUG_EX
   bool dont_store_a_structure = false;
   if (nlhs > 0)
     {
-      plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
-      pind = mxGetPr(plhs[0]);
-      if (no_error)
-        pind[0] = 0;
-      else
-        pind[0] = 1;
-      if (nlhs > 1)
+      if (block >= 0)
         {
-          if (block >= 0)
+          if (evaluate)
             {
-              if (evaluate)
-                {
-                  vector<double> residual = interprete.get_residual();
-                  plhs[1] = mxCreateDoubleMatrix(int (residual.size()/double (col_y)), int (col_y), mxREAL);
-                  pind = mxGetPr(plhs[1]);
-                  for (i = 0; i < residual.size(); i++)
-                    pind[i] = residual[i];
-                }
-              else
-                {
-                  int out_periods;
-                  if (extended_path)
-                    out_periods = max_periods + y_kmin;
-                  else
-                    out_periods = row_y;
-                  plhs[1] = mxCreateDoubleMatrix(out_periods, int (col_y), mxREAL);
-                  pind = mxGetPr(plhs[1]);
-                  for (i = 0; i < out_periods*col_y; i++)
-                    pind[i] = y[i];
-                }
+              vector<double> residual = interprete.get_residual();
+              plhs[0] = mxCreateDoubleMatrix(int (residual.size()/double (col_y)), int (col_y), mxREAL);
+              pind = mxGetPr(plhs[0]);
+              for (i = 0; i < residual.size(); i++)
+                pind[i] = residual[i];
             }
           else
             {
@@ -1159,99 +1106,111 @@ main(int nrhs, const char *prhs[])
               if (extended_path)
                 out_periods = max_periods + y_kmin;
               else
-                out_periods = col_y;
-              plhs[1] = mxCreateDoubleMatrix(int (row_y), out_periods, mxREAL);
-              pind = mxGetPr(plhs[1]);
-              if (evaluate)
+                out_periods = row_y;
+              plhs[0] = mxCreateDoubleMatrix(out_periods, int (col_y), mxREAL);
+              pind = mxGetPr(plhs[0]);
+              for (i = 0; i < out_periods*col_y; i++)
+                pind[i] = y[i];
+            }
+        }
+      else
+        {
+          int out_periods;
+          if (extended_path)
+            out_periods = max_periods + y_kmin;
+          else
+            out_periods = col_y;
+          plhs[0] = mxCreateDoubleMatrix(int (row_y), out_periods, mxREAL);
+          pind = mxGetPr(plhs[0]);
+          if (evaluate)
+            {
+              vector<double> residual = interprete.get_residual();
+              for (i = 0; i < residual.size(); i++)
+                pind[i] = residual[i];
+            }
+          else
+            for (i = 0; i < row_y*out_periods; i++)
+              pind[i] = y[i];
+        }
+      if (nlhs > 1)
+        {
+          if (evaluate)
+            {
+              int jacob_field_number = 0, jacob_exo_field_number = 0, jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0;
+              if (!block_structur)
+                {
+                  const char *field_names[] = {"g1", "g1_x", "g1_xd", "g1_o"};
+                  jacob_field_number = 0;
+                  jacob_exo_field_number = 1;
+                  jacob_exo_det_field_number = 2;
+                  jacob_other_endo_field_number = 3;
+                  mwSize dims[1] = {static_cast<mwSize>(nb_blocks) };
+                  plhs[1] = mxCreateStructArray(1, dims, 4, field_names);
+                }
+              else if (!mxIsStruct(block_structur))
                 {
-                  vector<double> residual = interprete.get_residual();
-                  for (i = 0; i < residual.size(); i++)
-                    pind[i] = residual[i];
+                  plhs[1] = interprete.get_jacob(0);
+                  //mexCallMATLAB(0,NULL, 1, &plhs[1], "disp");
+                  dont_store_a_structure = true;
                 }
               else
-                for (i = 0; i < row_y*out_periods; i++)
-                  pind[i] = y[i];
-            }
-          if (nlhs > 2)
-            {
-              if (evaluate)
                 {
-                  int jacob_field_number = 0, jacob_exo_field_number = 0, jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0;
-                  if (!block_structur)
-                    {
-                      const char *field_names[] = {"g1", "g1_x", "g1_xd", "g1_o"};
-                      jacob_field_number = 0;
-                      jacob_exo_field_number = 1;
-                      jacob_exo_det_field_number = 2;
-                      jacob_other_endo_field_number = 3;
-                      mwSize dims[1] = {static_cast<mwSize>(nb_blocks) };
-                      plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
-                    }
-                  else if (!mxIsStruct(block_structur))
-                    {
-                      plhs[2] = interprete.get_jacob(0);
-                      //mexCallMATLAB(0,NULL, 1, &plhs[2], "disp");
-                      dont_store_a_structure = true;
-                    }
-                  else
-                    {
-                      plhs[2] = block_structur;
-                      jacob_field_number = mxAddField(plhs[2], "g1");
-                      if (jacob_field_number == -1)
-                        DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob to the structArray\n");
-                      jacob_exo_field_number = mxAddField(plhs[2], "g1_x");
-                      if (jacob_exo_field_number == -1)
-                        DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob_exo to the structArray\n");
-                      jacob_exo_det_field_number = mxAddField(plhs[2], "g1_xd");
-                      if (jacob_exo_det_field_number == -1)
-                        DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob_exo_det to the structArray\n");
-                      jacob_other_endo_field_number = mxAddField(plhs[2], "g1_o");
-                      if (jacob_other_endo_field_number == -1)
-                        DYN_MEX_FUNC_ERR_MSG_TXT("Fatal error in bytecode: in main, cannot add extra field jacob_other_endo to the structArray\n");
-                    }
-                  if (!dont_store_a_structure)
+                  plhs[1] = block_structur;
+                  jacob_field_number = mxAddField(plhs[1], "g1");
+                  if (jacob_field_number == -1)
+                    mexErrMsgTxt("Fatal error in bytecode: in main, cannot add extra field jacob to the structArray\n");
+                  jacob_exo_field_number = mxAddField(plhs[1], "g1_x");
+                  if (jacob_exo_field_number == -1)
+                    mexErrMsgTxt("Fatal error in bytecode: in main, cannot add extra field jacob_exo to the structArray\n");
+                  jacob_exo_det_field_number = mxAddField(plhs[1], "g1_xd");
+                  if (jacob_exo_det_field_number == -1)
+                    mexErrMsgTxt("Fatal error in bytecode: in main, cannot add extra field jacob_exo_det to the structArray\n");
+                  jacob_other_endo_field_number = mxAddField(plhs[1], "g1_o");
+                  if (jacob_other_endo_field_number == -1)
+                    mexErrMsgTxt("Fatal error in bytecode: in main, cannot add extra field jacob_other_endo to the structArray\n");
+                }
+              if (!dont_store_a_structure)
+                {
+                  for (int i = 0; i < nb_blocks; i++)
                     {
-                      for (int i = 0; i < nb_blocks; i++)
+                      mxSetFieldByNumber(plhs[1], i, jacob_field_number, interprete.get_jacob(i));
+                      if (!steady_state)
                         {
-                          mxSetFieldByNumber(plhs[2], i, jacob_field_number, interprete.get_jacob(i));
-                          if (!steady_state)
-                            {
-                              mxSetFieldByNumber(plhs[2], i, jacob_exo_field_number, interprete.get_jacob_exo(i));
-                              mxSetFieldByNumber(plhs[2], i, jacob_exo_det_field_number, interprete.get_jacob_exo_det(i));
-                              mxSetFieldByNumber(plhs[2], i, jacob_other_endo_field_number, interprete.get_jacob_other_endo(i));
-                            }
+                          mxSetFieldByNumber(plhs[1], i, jacob_exo_field_number, interprete.get_jacob_exo(i));
+                          mxSetFieldByNumber(plhs[1], i, jacob_exo_det_field_number, interprete.get_jacob_exo_det(i));
+                          mxSetFieldByNumber(plhs[1], i, jacob_other_endo_field_number, interprete.get_jacob_other_endo(i));
                         }
                     }
                 }
-              else
+            }
+          else
+            {
+              plhs[1] = mxCreateDoubleMatrix(int (row_x), int (col_x), mxREAL);
+              pind = mxGetPr(plhs[1]);
+              for (i = 0; i < row_x*col_x; i++)
                 {
-                  plhs[2] = mxCreateDoubleMatrix(int (row_x), int (col_x), mxREAL);
-                  pind = mxGetPr(plhs[2]);
-                  for (i = 0; i < row_x*col_x; i++)
-                    {
-                      pind[i] = x[i];
-                    }
-
+                  pind[i] = x[i];
                 }
+
+            }
+          if (nlhs > 2)
+            {
+              plhs[2] = mxCreateDoubleMatrix(int (row_y), int (col_y), mxREAL);
+              pind = mxGetPr(plhs[2]);
+              for (i = 0; i < row_y*col_y; i++)
+                pind[i] = y[i];
               if (nlhs > 3)
                 {
-                  plhs[3] = mxCreateDoubleMatrix(int (row_y), int (col_y), mxREAL);
+                  mxArray *GlobalTemporaryTerms = interprete.get_Temporary_Terms();
+                  size_t nb_temp_terms = mxGetM(GlobalTemporaryTerms);
+                  plhs[3] = mxCreateDoubleMatrix(int (nb_temp_terms), 1, mxREAL);
                   pind = mxGetPr(plhs[3]);
-                  for (i = 0; i < row_y*col_y; i++)
-                    pind[i] = y[i];
-                  if (nlhs > 4)
-                    {
-                      mxArray *GlobalTemporaryTerms = interprete.get_Temporary_Terms();
-                      size_t nb_temp_terms = mxGetM(GlobalTemporaryTerms);
-                      plhs[4] = mxCreateDoubleMatrix(int (nb_temp_terms), 1, mxREAL);
-                      pind = mxGetPr(plhs[4]);
-                      double *tt = mxGetPr(GlobalTemporaryTerms);
-                      for (i = 0; i < nb_temp_terms; i++)
-                        pind[i] = tt[i];
-                    }
+                  double *tt = mxGetPr(GlobalTemporaryTerms);
+                  for (i = 0; i < nb_temp_terms; i++)
+                    pind[i] = tt[i];
                 }
-
             }
+
         }
     }
 #else
diff --git a/mex/sources/bytecode/testing/mex_interface.hh b/mex/sources/bytecode/testing/mex_interface.hh
index 9e8989189341b8c63fb6a97ff7e1f45e0ac4ed8a..b8b0dbfb0ad7b82a07ec747423581792f82809fe 100644
--- a/mex/sources/bytecode/testing/mex_interface.hh
+++ b/mex/sources/bytecode/testing/mex_interface.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2017 Dynare Team
+ * Copyright © 2007-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -29,13 +29,6 @@
 #include <algorithm>
 using namespace std;
 
-#if !defined(DYN_MEX_FUNC_ERR_MSG_TXT)
-# define DYN_MEX_FUNC_ERR_MSG_TXT(str)          \
-  do {                                          \
-    mexPrintf("%s\n", str);                     \
-  } while (0)
-#endif
-
 typedef unsigned int mwIndex;
 typedef unsigned int mwSize;
 
diff --git a/mex/sources/dynare_simul_/dynare_simul_.cc b/mex/sources/dynare_simul_/dynare_simul_.cc
index 58a6396092a80d5673d015424093be14bc29cc8e..faaf146b0383e743b3441757547173a3ba33adbc 100644
--- a/mex/sources/dynare_simul_/dynare_simul_.cc
+++ b/mex/sources/dynare_simul_/dynare_simul_.cc
@@ -1,6 +1,6 @@
 /*
  * Copyright © 2005-2011 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -51,8 +51,8 @@ extern "C" {
   mexFunction(int nlhs, mxArray *plhs[],
               int nhrs, const mxArray *prhs[])
   {
-    if (nhrs != 12 || nlhs != 2)
-      DYN_MEX_FUNC_ERR_MSG_TXT("dynare_simul_ must have at exactly 12 input parameters and 2 output arguments.\n");
+    if (nhrs != 12 || nlhs != 1)
+      mexErrMsgTxt("dynare_simul_ must have at exactly 12 input parameters and 1 output argument.");
 
     int order = static_cast<int>(mxGetScalar(prhs[0]));
     int nstat = static_cast<int>(mxGetScalar(prhs[1]));
@@ -74,22 +74,22 @@ extern "C" {
 
     int ny = nstat + npred + nboth + nforw;
     if (ny != static_cast<int>(ystart_dim[0]))
-      DYN_MEX_FUNC_ERR_MSG_TXT("ystart has wrong number of rows.\n");
+      mexErrMsgTxt("ystart has wrong number of rows.\n");
     if (1 != ystart_dim[1])
-      DYN_MEX_FUNC_ERR_MSG_TXT("ystart has wrong number of cols.\n");
+      mexErrMsgTxt("ystart has wrong number of cols.\n");
     int nper = shocks_dim[1];
     if (nexog != static_cast<int>(shocks_dim[0]))
-      DYN_MEX_FUNC_ERR_MSG_TXT("shocks has a wrong number of rows.\n");
+      mexErrMsgTxt("shocks has a wrong number of rows.\n");
     if (nexog != static_cast<int>(vcov_dim[0]))
-      DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of rows.\n");
+      mexErrMsgTxt("vcov has a wrong number of rows.\n");
     if (nexog != static_cast<int>(vcov_dim[1]))
-      DYN_MEX_FUNC_ERR_MSG_TXT("vcov has a wrong number of cols.\n");
+      mexErrMsgTxt("vcov has a wrong number of cols.\n");
     if (ny != static_cast<int>(ysteady_dim[0]))
-      DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of rows.\n");
+      mexErrMsgTxt("ysteady has wrong number of rows.\n");
     if (1 != ysteady_dim[1])
-      DYN_MEX_FUNC_ERR_MSG_TXT("ysteady has wrong number of cols.\n");
+      mexErrMsgTxt("ysteady has wrong number of cols.\n");
 
-    plhs[1] = mxCreateDoubleMatrix(ny, nper, mxREAL);
+    plhs[0] = mxCreateDoubleMatrix(ny, nper, mxREAL);
 
     try
       {
@@ -102,13 +102,13 @@ extern "C" {
           {
             const mxArray *gk_m = mxGetField(dr, 0, ("g_" + std::to_string(dim)).c_str());
             if (!gk_m)
-              DYN_MEX_FUNC_ERR_MSG_TXT(("Can't find field `g_" + std::to_string(dim) + "' in structured passed as last argument").c_str());
+              mexErrMsgTxt(("Can't find field `g_" + std::to_string(dim) + "' in structured passed as last argument").c_str());
             ConstTwoDMatrix gk(gk_m);
             FFSTensor ft(ny, npred+nboth+nexog, dim);
             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());
+              mexErrMsgTxt(("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());
+              mexErrMsgTxt(("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();
             ft.add(1.0, gk);
             pol.insert(std::make_unique<UFSTensor>(ft));
@@ -123,21 +123,20 @@ extern "C" {
         // simulate and copy the results
         TwoDMatrix res_mat{dr.simulate(DecisionRule::emethod::horner, nper,
                                        ConstVector{ystart}, sr)};
-        TwoDMatrix res_tmp_mat{plhs[1]};
+        TwoDMatrix res_tmp_mat{plhs[0]};
         res_tmp_mat = const_cast<const TwoDMatrix &>(res_mat);
       }
     catch (const KordException &e)
       {
-        DYN_MEX_FUNC_ERR_MSG_TXT("Caught Kord exception.");
+        mexErrMsgTxt("Caught Kord exception.");
       }
     catch (const TLException &e)
       {
-        DYN_MEX_FUNC_ERR_MSG_TXT("Caught TL exception.");
+        mexErrMsgTxt("Caught TL exception.");
       }
     catch (SylvException &e)
       {
-        DYN_MEX_FUNC_ERR_MSG_TXT("Caught Sylv exception.");
+        mexErrMsgTxt("Caught Sylv exception.");
       }
-    plhs[0] = mxCreateDoubleScalar(0);
   }
 };
diff --git a/mex/sources/dynmex.h b/mex/sources/dynmex.h
index 52a71bef5a9e3df359b726d3f791906c908d3a10..4bb0ca7391278d0b7f38bdbc65880b24ec251750 100644
--- a/mex/sources/dynmex.h
+++ b/mex/sources/dynmex.h
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2009-2019 Dynare Team
+ * Copyright © 2009-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -26,20 +26,6 @@
 
 #include <mex.h>
 
-/*
- * Fix for trac ticket Ticket #137
- */
-#if !defined(DYN_MEX_FUNC_ERR_MSG_TXT)
-# define DYN_MEX_FUNC_ERR_MSG_TXT(str)          \
-  do {                                          \
-    mexPrintf("%s\n", str);                     \
-    int i;                                      \
-    for (i = 0; i < nlhs; i++)                  \
-      plhs[i] = mxCreateDoubleScalar(1);        \
-    return;                                     \
-  } while (0)
-#endif
-
 #if defined(MATLAB_MEX_FILE) && MATLAB_VERSION < 0x0805
 # define mxIsScalar(x) (mxGetM(x) == 1 && mxGetN(x) == 1)
 #endif
diff --git a/mex/sources/gensylv/gensylv.cc b/mex/sources/gensylv/gensylv.cc
index 15a980ce7b11d626ed3627303be85704cd6df659..c71f295f92fd3c4a0dc20c10154140f76344f226 100644
--- a/mex/sources/gensylv/gensylv.cc
+++ b/mex/sources/gensylv/gensylv.cc
@@ -1,6 +1,6 @@
 /*
  * Copyright © 2005-2011 Ondra Kamenik
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -51,8 +51,8 @@ extern "C" {
   mexFunction(int nlhs, mxArray *plhs[],
               int nhrs, const mxArray *prhs[])
   {
-    if (nhrs != 5 || nlhs > 3 || nlhs < 2)
-      DYN_MEX_FUNC_ERR_MSG_TXT("Gensylv: Must have exactly 5 input args and either 2 or 3 output args.");
+    if (nhrs != 5 || nlhs > 2 || nlhs < 1)
+      mexErrMsgTxt("Gensylv: Must have exactly 5 input args and either 1 or 2 output args.");
 
     auto order = static_cast<int>(mxGetScalar(prhs[0]));
     const mxArray *const A = prhs[1];
@@ -65,17 +65,17 @@ extern "C" {
     const mwSize *const Ddims = mxGetDimensions(D);
 
     if (Adims[0] != Adims[1])
-      DYN_MEX_FUNC_ERR_MSG_TXT("Matrix A must be a square matrix.");
+      mexErrMsgTxt("Matrix A must be a square matrix.");
     if (Adims[0] != Bdims[0])
-      DYN_MEX_FUNC_ERR_MSG_TXT("Matrix A and matrix B must have the same number of rows.");
+      mexErrMsgTxt("Matrix A and matrix B must have the same number of rows.");
     if (Adims[0] != Ddims[0])
-      DYN_MEX_FUNC_ERR_MSG_TXT("Matrix A and matrix B must have the same number of rows.");
+      mexErrMsgTxt("Matrix A and matrix B must have the same number of rows.");
     if (Cdims[0] != Cdims[1])
-      DYN_MEX_FUNC_ERR_MSG_TXT("Matrix C must be square.");
+      mexErrMsgTxt("Matrix C must be square.");
     if (Bdims[0] < Bdims[1])
-      DYN_MEX_FUNC_ERR_MSG_TXT("Matrix B must not have more columns than rows.");
+      mexErrMsgTxt("Matrix B must not have more columns than rows.");
     if (Ddims[1] != static_cast<mwSize>(power(Cdims[0], order)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("Matrix D has wrong number of columns.");
+      mexErrMsgTxt("Matrix D has wrong number of columns.");
 
     auto n = static_cast<int>(Adims[0]);
     auto m = static_cast<int>(Cdims[0]);
@@ -88,16 +88,15 @@ extern "C" {
     // solve (or solve and check)
     try
       {
-        if (nlhs == 2)
+        if (nlhs == 1)
           gen_sylv_solve(order, n, m, zero_cols, Avec, Bvec, Cvec, Xvec);
-        else if (nlhs == 3)
-          plhs[2] = gen_sylv_solve_and_check(order, n, m, zero_cols, Avec, Bvec, Cvec, Dvec, Xvec);
+        else if (nlhs == 2)
+          plhs[1] = gen_sylv_solve_and_check(order, n, m, zero_cols, Avec, Bvec, Cvec, Dvec, Xvec);
       }
     catch (const SylvException &e)
       {
-        DYN_MEX_FUNC_ERR_MSG_TXT(e.getMessage().c_str());
+        mexErrMsgTxt(e.getMessage().c_str());
       }
-    plhs[1] = X;
-    plhs[0] = mxCreateDoubleScalar(0);
+    plhs[0] = X;
   }
 };
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 0a44e820623286ba85237b3fe81be093d81ae3fc..4ae1151db72937c7deccb497838701d172c7a75b 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2008-2019 Dynare Team
+ * Copyright © 2008-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -82,8 +82,8 @@ extern "C" {
   mexFunction(int nlhs, mxArray *plhs[],
               int nrhs, const mxArray *prhs[])
   {
-    if (nrhs < 3 || nlhs < 2 || nlhs > 3)
-      DYN_MEX_FUNC_ERR_MSG_TXT("Must have at least 3 input parameters and takes 2 or 3 output parameters.");
+    if (nrhs < 3 || nlhs < 1 || nlhs > 2)
+      mexErrMsgTxt("Must have at least 3 input parameters and takes 1 or 2 output parameters.");
 
     // Give explicit names to input arguments
     const mxArray *dr_mx = prhs[0];
@@ -101,11 +101,11 @@ extern "C" {
     // Extract various fields from options_
     const int kOrder = get_int_field(options_mx, "order");
     if (kOrder < 1)
-      DYN_MEX_FUNC_ERR_MSG_TXT("options_.order must be at least 1");
+      mexErrMsgTxt("options_.order must be at least 1");
 
     const mxArray *use_dll_mx = mxGetField(options_mx, 0, "use_dll");
     if (!(use_dll_mx && mxIsLogicalScalar(use_dll_mx)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("options_.use_dll should be a logical scalar");
+      mexErrMsgTxt("options_.use_dll should be a logical scalar");
     bool use_dll = static_cast<bool>(mxGetScalar(use_dll_mx));
 
     double qz_criterium = 1+1e-6;
@@ -115,31 +115,31 @@ extern "C" {
 
     const mxArray *threads_mx = mxGetField(options_mx, 0, "threads");
     if (!threads_mx)
-      DYN_MEX_FUNC_ERR_MSG_TXT("Can't find field options_.threads");
+      mexErrMsgTxt("Can't find field options_.threads");
     const mxArray *num_threads_mx = mxGetField(threads_mx, 0, "k_order_perturbation");
     if (!(num_threads_mx && mxIsScalar(num_threads_mx) && mxIsNumeric(num_threads_mx)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("options_.threads.k_order_perturbation be a numeric scalar");
+      mexErrMsgTxt("options_.threads.k_order_perturbation be a numeric scalar");
     int num_threads = static_cast<int>(mxGetScalar(num_threads_mx));
 
     // Extract various fields from M_
     const mxArray *fname_mx = mxGetField(M_mx, 0, "fname");
     if (!(fname_mx && mxIsChar(fname_mx) && mxGetM(fname_mx) == 1))
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.fname should be a character string");
+      mexErrMsgTxt("M_.fname should be a character string");
     std::string fname{mxArrayToString(fname_mx)};
 
     const mxArray *params_mx = mxGetField(M_mx, 0, "params");
     if (!(params_mx && mxIsDouble(params_mx)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.params should be a double precision array");
+      mexErrMsgTxt("M_.params should be a double precision array");
     Vector modParams{ConstVector{params_mx}};
     if (!modParams.isFinite())
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.params contains NaN or Inf");
+      mexErrMsgTxt("M_.params contains NaN or Inf");
 
     const mxArray *sigma_e_mx = mxGetField(M_mx, 0, "Sigma_e");
     if (!(sigma_e_mx && mxIsDouble(sigma_e_mx) && mxGetM(sigma_e_mx) == mxGetN(sigma_e_mx)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.Sigma_e should be a double precision square matrix");
+      mexErrMsgTxt("M_.Sigma_e should be a double precision square matrix");
     TwoDMatrix vCov{ConstTwoDMatrix{sigma_e_mx}};
     if (!vCov.isFinite())
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.Sigma_e contains NaN or Inf");
+      mexErrMsgTxt("M_.Sigma_e contains NaN or Inf");
 
     const int nStat = get_int_field(M_mx, "nstatic");
     const int nPred = get_int_field(M_mx, "npred");
@@ -153,42 +153,42 @@ extern "C" {
     const mxArray *lead_lag_incidence_mx = mxGetField(M_mx, 0, "lead_lag_incidence");
     if (!(lead_lag_incidence_mx && mxIsDouble(lead_lag_incidence_mx) && mxGetM(lead_lag_incidence_mx) == 3
           && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(nEndo)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.lead_lag_incidence should be a double precision matrix with 3 rows and M_.endo_nbr columns");
+      mexErrMsgTxt("M_.lead_lag_incidence should be a double precision matrix with 3 rows and M_.endo_nbr columns");
     ConstTwoDMatrix llincidence{lead_lag_incidence_mx};
 
     const mxArray *nnzderivatives_mx = mxGetField(M_mx, 0, "NNZDerivatives");
     if (!(nnzderivatives_mx && mxIsDouble(nnzderivatives_mx)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.NNZDerivatives should be a double precision array");
+      mexErrMsgTxt("M_.NNZDerivatives should be a double precision array");
     ConstVector NNZD{nnzderivatives_mx};
     if (NNZD.length() < kOrder || NNZD[kOrder-1] == -1)
-      DYN_MEX_FUNC_ERR_MSG_TXT("The derivatives were not computed for the required order. Make sure that you used the right order option inside the `stoch_simul' command");
+      mexErrMsgTxt("The derivatives were not computed for the required order. Make sure that you used the right order option inside the `stoch_simul' command");
 
     const mxArray *endo_names_mx = mxGetField(M_mx, 0, "endo_names");
     if (!(endo_names_mx && mxIsCell(endo_names_mx) && mxGetNumberOfElements(endo_names_mx) == static_cast<size_t>(nEndo)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.endo_names should be a cell array of M_.endo_nbr elements");
+      mexErrMsgTxt("M_.endo_names should be a cell array of M_.endo_nbr elements");
     std::vector<std::string> endoNames = DynareMxArrayToString(endo_names_mx);
 
     const mxArray *exo_names_mx = mxGetField(M_mx, 0, "exo_names");
     if (!(exo_names_mx && mxIsCell(exo_names_mx) && mxGetNumberOfElements(exo_names_mx) == static_cast<size_t>(nExog)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.exo_names should be a cell array of M_.exo_nbr elements");
+      mexErrMsgTxt("M_.exo_names should be a cell array of M_.exo_nbr elements");
     std::vector<std::string> exoNames = DynareMxArrayToString(exo_names_mx);
 
     const mxArray *dynamic_tmp_nbr_mx = mxGetField(M_mx, 0, "dynamic_tmp_nbr");
     if (!(dynamic_tmp_nbr_mx && mxIsDouble(dynamic_tmp_nbr_mx) && mxGetNumberOfElements(dynamic_tmp_nbr_mx) >= static_cast<size_t>(kOrder+1)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("M_.dynamic_tmp_nbr should be a double precision array with strictly more elements than the order of derivation");
+      mexErrMsgTxt("M_.dynamic_tmp_nbr should be a double precision array with strictly more elements than the order of derivation");
     int ntt = std::accumulate(mxGetPr(dynamic_tmp_nbr_mx), mxGetPr(dynamic_tmp_nbr_mx)+kOrder+1, 0);
 
     // Extract various fields from dr
     const mxArray *ys_mx = mxGetField(dr_mx, 0, "ys"); // and not in order of dr.order_var
     if (!(ys_mx && mxIsDouble(ys_mx)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("dr.ys should be a double precision array");
+      mexErrMsgTxt("dr.ys should be a double precision array");
     Vector ySteady{ConstVector{ys_mx}};
     if (!ySteady.isFinite())
-      DYN_MEX_FUNC_ERR_MSG_TXT("dr.ys contains NaN or Inf");
+      mexErrMsgTxt("dr.ys contains NaN or Inf");
 
     const mxArray *order_var_mx = mxGetField(dr_mx, 0, "order_var");
     if (!(order_var_mx && mxIsDouble(order_var_mx) && mxGetNumberOfElements(order_var_mx) == static_cast<size_t>(nEndo)))
-      DYN_MEX_FUNC_ERR_MSG_TXT("dr.order_var should be a double precision array of M_.endo_nbr elements");
+      mexErrMsgTxt("dr.order_var should be a double precision array of M_.endo_nbr elements");
     std::vector<int> dr_order(nEndo);
     std::transform(mxGetPr(order_var_mx), mxGetPr(order_var_mx)+nEndo, dr_order.begin(),
                    [](double x) { return static_cast<int>(x)-1; });
@@ -241,7 +241,7 @@ extern "C" {
         const char *g_fieldnames_c[kOrder+1];
         for (int i = 0; i <= kOrder; i++)
           g_fieldnames_c[i] = g_fieldnames[i].c_str();
-        plhs[1] = mxCreateStructMatrix(1, 1, kOrder+1, g_fieldnames_c);
+        plhs[0] = mxCreateStructMatrix(1, 1, kOrder+1, g_fieldnames_c);
 
         // Fill that structure
         for (int i = 0; i <= kOrder; i++)
@@ -251,10 +251,10 @@ extern "C" {
             const ConstVector &vec = t.getData();
             assert(vec.skip() == 1);
             std::copy_n(vec.base(), vec.length(), mxGetPr(tmp));
-            mxSetField(plhs[1], 0, ("g_" + std::to_string(i)).c_str(), tmp);
+            mxSetField(plhs[0], 0, ("g_" + std::to_string(i)).c_str(), tmp);
           }
 
-        if (nlhs > 2)
+        if (nlhs > 1)
           {
             /* Return as 3rd argument a struct containing derivatives in Dynare
                format (unfolded matrices, without Taylor coefficient) up to 3rd
@@ -264,51 +264,50 @@ extern "C" {
             size_t nfields = (kOrder == 1 ? 2 : (kOrder == 2 ? 6 : 12));
             const char *c_fieldnames[] = { "gy", "gu", "gyy", "gyu", "guu", "gss",
                                            "gyyy", "gyyu", "gyuu", "guuu", "gyss", "guss" };
-            plhs[2] = mxCreateStructMatrix(1, 1, nfields, c_fieldnames);
+            plhs[1] = mxCreateStructMatrix(1, 1, nfields, c_fieldnames);
 
-            copy_derivatives(plhs[2], Symmetry{1, 0, 0, 0}, derivs, "gy");
-            copy_derivatives(plhs[2], Symmetry{0, 1, 0, 0}, derivs, "gu");
+            copy_derivatives(plhs[1], Symmetry{1, 0, 0, 0}, derivs, "gy");
+            copy_derivatives(plhs[1], Symmetry{0, 1, 0, 0}, derivs, "gu");
             if (kOrder >= 2)
               {
-                copy_derivatives(plhs[2], Symmetry{2, 0, 0, 0}, derivs, "gyy");
-                copy_derivatives(plhs[2], Symmetry{0, 2, 0, 0}, derivs, "guu");
-                copy_derivatives(plhs[2], Symmetry{1, 1, 0, 0}, derivs, "gyu");
-                copy_derivatives(plhs[2], Symmetry{0, 0, 0, 2}, derivs, "gss");
+                copy_derivatives(plhs[1], Symmetry{2, 0, 0, 0}, derivs, "gyy");
+                copy_derivatives(plhs[1], Symmetry{0, 2, 0, 0}, derivs, "guu");
+                copy_derivatives(plhs[1], Symmetry{1, 1, 0, 0}, derivs, "gyu");
+                copy_derivatives(plhs[1], Symmetry{0, 0, 0, 2}, derivs, "gss");
               }
             if (kOrder >= 3)
               {
-                copy_derivatives(plhs[2], Symmetry{3, 0, 0, 0}, derivs, "gyyy");
-                copy_derivatives(plhs[2], Symmetry{0, 3, 0, 0}, derivs, "guuu");
-                copy_derivatives(plhs[2], Symmetry{2, 1, 0, 0}, derivs, "gyyu");
-                copy_derivatives(plhs[2], Symmetry{1, 2, 0, 0}, derivs, "gyuu");
-                copy_derivatives(plhs[2], Symmetry{1, 0, 0, 2}, derivs, "gyss");
-                copy_derivatives(plhs[2], Symmetry{0, 1, 0, 2}, derivs, "guss");
+                copy_derivatives(plhs[1], Symmetry{3, 0, 0, 0}, derivs, "gyyy");
+                copy_derivatives(plhs[1], Symmetry{0, 3, 0, 0}, derivs, "guuu");
+                copy_derivatives(plhs[1], Symmetry{2, 1, 0, 0}, derivs, "gyyu");
+                copy_derivatives(plhs[1], Symmetry{1, 2, 0, 0}, derivs, "gyuu");
+                copy_derivatives(plhs[1], Symmetry{1, 0, 0, 2}, derivs, "gyss");
+                copy_derivatives(plhs[1], Symmetry{0, 1, 0, 2}, derivs, "guss");
               }
           }
       }
     catch (const KordException &e)
       {
         e.print();
-        DYN_MEX_FUNC_ERR_MSG_TXT(("dynare:k_order_perturbation: Caught Kord exception: " + e.get_message()).c_str());
+        mexErrMsgTxt(("dynare:k_order_perturbation: Caught Kord exception: " + e.get_message()).c_str());
       }
     catch (const TLException &e)
       {
         e.print();
-        DYN_MEX_FUNC_ERR_MSG_TXT("dynare:k_order_perturbation: Caught TL exception");
+        mexErrMsgTxt("dynare:k_order_perturbation: Caught TL exception");
       }
     catch (SylvException &e)
       {
         e.printMessage();
-        DYN_MEX_FUNC_ERR_MSG_TXT("dynare:k_order_perturbation: Caught Sylv exception");
+        mexErrMsgTxt("dynare:k_order_perturbation: Caught Sylv exception");
       }
     catch (const DynareException &e)
       {
-        DYN_MEX_FUNC_ERR_MSG_TXT(("dynare:k_order_perturbation: Caught KordDynare exception: " + e.message()).c_str());
+        mexErrMsgTxt(("dynare:k_order_perturbation: Caught KordDynare exception: " + e.message()).c_str());
       }
     catch (const ogu::Exception &e)
       {
-        DYN_MEX_FUNC_ERR_MSG_TXT(("dynare:k_order_perturbation: Caught general exception: " + e.message()).c_str());
+        mexErrMsgTxt(("dynare:k_order_perturbation: Caught general exception: " + e.message()).c_str());
       }
-    plhs[0] = mxCreateDoubleScalar(0);
   } // end of mexFunction()
 } // end of extern C
diff --git a/mex/sources/kalman_steady_state/kalman_steady_state.cc b/mex/sources/kalman_steady_state/kalman_steady_state.cc
index 0a0db3181882d52b5eed5652ba451d8713257e4c..14b61015b34e4ead9bb980c70cb30f296af684a6 100644
--- a/mex/sources/kalman_steady_state/kalman_steady_state.cc
+++ b/mex/sources/kalman_steady_state/kalman_steady_state.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2009-2019 Dynare Team.
+ * Copyright © 2009-2020 Dynare Team.
  *
  * This file is part of Dynare.
  *
@@ -76,10 +76,10 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   // Check the number of arguments and set some flags.
   bool measurement_error_flag = true;
   if (nrhs < 3 || 4 < nrhs)
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state accepts either 3 or 4 input arguments!");
+    mexErrMsgTxt("kalman_steady_state accepts either 3 or 4 input arguments!");
 
-  if (nlhs < 1 || 2 < nlhs)
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state requires at least 1, but no more than 2, output arguments!");
+  if (nlhs != 1)
+    mexErrMsgTxt("kalman_steady_state accepts exactly one output argument!");
 
   if (nrhs == 3)
     measurement_error_flag = false;
@@ -87,30 +87,30 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   // Check the type of the input arguments and get the size of the matrices.
   lapack_int n = mxGetM(prhs[0]);
   if (static_cast<size_t>(n) != mxGetN(prhs[0]))
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The first input argument (T) must be a square matrix!");
+    mexErrMsgTxt("kalman_steady_state: The first input argument (T) must be a square matrix!");
   if (mxIsNumeric(prhs[0]) == 0 || mxIsComplex(prhs[0]) == 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The first input argument (T) must be a real matrix!");
+    mexErrMsgTxt("kalman_steady_state: The first input argument (T) must be a real matrix!");
 
   lapack_int q = mxGetM(prhs[1]);
   if (static_cast<size_t>(q) != mxGetN(prhs[1]))
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The second input argument (QQ) must be a square matrix!");
+    mexErrMsgTxt("kalman_steady_state: The second input argument (QQ) must be a square matrix!");
   if (mxIsNumeric(prhs[1]) == 0 || mxIsComplex(prhs[1]) == 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The second input argument (QQ) must be a real matrix!");
+    mexErrMsgTxt("kalman_steady_state: The second input argument (QQ) must be a real matrix!");
   if (q != n)
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The size of the second input argument (QQ) must match the size of the first argument (T)!");
+    mexErrMsgTxt("kalman_steady_state: The size of the second input argument (QQ) must match the size of the first argument (T)!");
   lapack_int p = mxGetN(prhs[2]);
   if (mxGetM(prhs[2]) != static_cast<size_t>(n))
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The number of rows of the third argument (Z) must match the number of rows of the first argument (T)!");
+    mexErrMsgTxt("kalman_steady_state: The number of rows of the third argument (Z) must match the number of rows of the first argument (T)!");
   if (mxIsNumeric(prhs[2]) == 0 || mxIsComplex(prhs[2]) == 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The third input argument (Z) must be a real matrix!");
+    mexErrMsgTxt("kalman_steady_state: The third input argument (Z) must be a real matrix!");
   if (measurement_error_flag)
     {
       if (mxGetM(prhs[3]) != mxGetN(prhs[3]))
-        DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The fourth input argument (H) must be a square matrix!");
+        mexErrMsgTxt("kalman_steady_state: The fourth input argument (H) must be a square matrix!");
       if (mxGetM(prhs[3]) != static_cast<size_t>(p))
-        DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The number of rows of the fourth input argument (H) must match the number of rows of the third input argument!");
+        mexErrMsgTxt("kalman_steady_state: The number of rows of the fourth input argument (H) must match the number of rows of the third input argument!");
       if (mxIsNumeric(prhs[3]) == 0 || mxIsComplex(prhs[3]) == 1)
-        DYN_MEX_FUNC_ERR_MSG_TXT("kalman_steady_state: The fifth input argument (H) must be a real matrix!");
+        mexErrMsgTxt("kalman_steady_state: The fifth input argument (H) must be a real matrix!");
     }
   // Get input matrices.
   const double *T = mxGetPr(prhs[0]);
@@ -147,8 +147,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   auto DWORK = std::make_unique<double[]>(LDWORK);
   auto BWORK = std::make_unique<lapack_int[]>(nn);
   // Initialize the output of the mex file
-  plhs[1] = mxCreateDoubleMatrix(n, n, mxREAL);
-  double *P = mxGetPr(plhs[1]);
+  plhs[0] = mxCreateDoubleMatrix(n, n, mxREAL);
+  double *P = mxGetPr(plhs[0]);
   // Call the slicot routine
   sb02od("D", // We want to solve a discrete Riccati equation.
          "B", // Matrices Z and H are given.
@@ -164,25 +164,24 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     case 0:
       break;
     case 1:
-      DYN_MEX_FUNC_ERR_MSG_TXT("The computed extended matrix pencil is singular, possibly due to rounding errors.\n");
+      mexErrMsgTxt("The computed extended matrix pencil is singular, possibly due to rounding errors.");
       break;
     case 2:
-      DYN_MEX_FUNC_ERR_MSG_TXT("The QZ (or QR) algorithm failed!\n");
+      mexErrMsgTxt("The QZ (or QR) algorithm failed!");
       break;
     case 3:
-      DYN_MEX_FUNC_ERR_MSG_TXT("The reordering of the (generalized) eigenvalues failed!\n");
+      mexErrMsgTxt("The reordering of the (generalized) eigenvalues failed!");
       break;
     case 4:
-      DYN_MEX_FUNC_ERR_MSG_TXT("After reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues\n in the (generalized) Schur form no longer satisfy the stability condition; this could also be caused due to scaling.");
+      mexErrMsgTxt("After reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues\n in the (generalized) Schur form no longer satisfy the stability condition; this could also be caused due to scaling.");
       break;
     case 5:
-      DYN_MEX_FUNC_ERR_MSG_TXT("The computed dimension of the solution does not equal n!\n");
+      mexErrMsgTxt("The computed dimension of the solution does not equal n!");
       break;
     case 6:
-      DYN_MEX_FUNC_ERR_MSG_TXT("A singular matrix was encountered during the computation of the solution matrix P!\n");
+      mexErrMsgTxt("A singular matrix was encountered during the computation of the solution matrix P!");
       break;
     default:
-      DYN_MEX_FUNC_ERR_MSG_TXT("Unknown problem!\n");
+      mexErrMsgTxt("Unknown problem!");
     }
-  plhs[0] = mxCreateDoubleScalar(0);
 }
diff --git a/mex/sources/kronecker/A_times_B_kronecker_C.cc b/mex/sources/kronecker/A_times_B_kronecker_C.cc
index 071122e3f9932975d72d3e50a7aec77a5fe5fca4..95d0f72680ebc370c33c5c53f207817ee6512388 100644
--- a/mex/sources/kronecker/A_times_B_kronecker_C.cc
+++ b/mex/sources/kronecker/A_times_B_kronecker_C.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2019 Dynare Team
+ * Copyright © 2007-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -68,8 +68,11 @@ void
 mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
   // Check input and output:
-  if (nrhs > 3 || nrhs < 2)
-    DYN_MEX_FUNC_ERR_MSG_TXT("A_times_B_kronecker_C takes 2 or 3 input arguments and provides 2 output arguments.");
+  if (nrhs > 3 || nrhs < 2 || nlhs != 1)
+    {
+      mexErrMsgTxt("A_times_B_kronecker_C takes 2 or 3 input arguments and provides 1 output argument.");
+      return; // Needed to shut up some GCC warnings
+    }
 
   // Get & Check dimensions (columns and rows):
   size_t mA = mxGetM(prhs[0]);
@@ -82,12 +85,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
       mC = mxGetM(prhs[2]);
       nC = mxGetN(prhs[2]);
       if (mB*mC != nA)
-        DYN_MEX_FUNC_ERR_MSG_TXT("Input dimension error!");
+        mexErrMsgTxt("Input dimension error!");
     }
   else // A·(B⊗B) is to be computed.
     {
       if (mB*mB != nA)
-        DYN_MEX_FUNC_ERR_MSG_TXT("Input dimension error!");
+        mexErrMsgTxt("Input dimension error!");
     }
   // Get input matrices:
   const double *A = mxGetPr(prhs[0]);
@@ -108,6 +111,4 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     full_A_times_kronecker_B_B(A, B, D, mA, nA, mB, nB);
   else
     full_A_times_kronecker_B_C(A, B, C, D, mA, nA, mB, nB, mC, nC);
-
-  plhs[1] = mxCreateDoubleScalar(0);
 }
diff --git a/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc b/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc
index 93265998cc32f06afd7a439677f4100a5fa13d81..5eaa7eee597f690ee86df9314b44348c0530617e 100644
--- a/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc
+++ b/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2007-2019 Dynare Team
+ * Copyright © 2007-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -135,11 +135,14 @@ void
 mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
   // Check input and output:
-  if (nrhs > 4 || nrhs < 3)
-    DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C takes 3 or 4 input arguments and provides 2 output arguments.");
+  if (nrhs > 4 || nrhs < 3 || nlhs != 1)
+    {
+      mexErrMsgTxt("sparse_hessian_times_B_kronecker_C takes 3 or 4 input arguments and provides 1 output argument.");
+      return; // Needed to shut up some GCC warnings
+    }
 
   if (!mxIsSparse(prhs[0]))
-    DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C: First input must be a sparse (dynare) hessian matrix.");
+    mexErrMsgTxt("sparse_hessian_times_B_kronecker_C: First input must be a sparse (dynare) hessian matrix.");
 
   // Get & Check dimensions (columns and rows):
   size_t mA = mxGetM(prhs[0]);
@@ -152,12 +155,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
       mC = mxGetM(prhs[2]);
       nC = mxGetN(prhs[2]);
       if (mB*mC != nA)
-        DYN_MEX_FUNC_ERR_MSG_TXT("Input dimension error!");
+        mexErrMsgTxt("Input dimension error!");
     }
   else // A·(B⊗B) is to be computed.
     {
       if (mB*mB != nA)
-        DYN_MEX_FUNC_ERR_MSG_TXT("Input dimension error!");
+        mexErrMsgTxt("Input dimension error!");
     }
   // Get input matrices:
   int numthreads;
@@ -185,6 +188,4 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     sparse_hessian_times_B_kronecker_B(isparseA, jsparseA, vsparseA, B, D, mA, nA, mB, nB, numthreads);
   else
     sparse_hessian_times_B_kronecker_C(isparseA, jsparseA, vsparseA, B, C, D, mA, nA, mB, nB, mC, nC, numthreads);
-
-  plhs[1] = mxCreateDoubleScalar(0);
 }
diff --git a/mex/sources/mjdgges/mjdgges.F08 b/mex/sources/mjdgges/mjdgges.F08
index dbc5da9d1246a4abdc3e9426b9633a4937da7ed0..5209c01d1d3340772f10e3cfeb063092aff307a7 100644
--- a/mex/sources/mjdgges/mjdgges.F08
+++ b/mex/sources/mjdgges/mjdgges.F08
@@ -60,8 +60,8 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   real(real64), dimension(:), allocatable, target :: vsl_target
 #endif
 
-  if (nrhs < 2 .or. nrhs > 4 .or. nlhs /= 7) then
-     call mexErrMsgTxt("MJDGGES: takes 2, 3 or 4 input arguments and exactly 7 output arguments.")
+  if (nrhs < 2 .or. nrhs > 4 .or. nlhs /= 6) then
+     call mexErrMsgTxt("MJDGGES: takes 2, 3 or 4 input arguments and exactly 6 output arguments.")
      return
   end if
 
@@ -95,24 +95,24 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
      zhreshold = 1e-6_real64
   end if
 
+  plhs(1) = mxCreateDoubleMatrix(n1, n1, mxREAL)
   plhs(2) = mxCreateDoubleMatrix(n1, n1, mxREAL)
   plhs(3) = mxCreateDoubleMatrix(n1, n1, mxREAL)
-  plhs(4) = mxCreateDoubleMatrix(n1, n1, mxREAL)
-  plhs(5) = mxCreateDoubleMatrix(1_mwSize, 1_mwSize, mxREAL)
-  plhs(6) = mxCreateDoubleMatrix(n1, 1_mwSize, mxCOMPLEX)
-  plhs(7) = mxCreateDoubleMatrix(1_mwSize, 1_mwSize, mxREAL)
-
-  s => mxGetPr(plhs(2))
-  t => mxGetPr(plhs(3))
-  sdim => mxGetPr(plhs(5))
+  plhs(4) = mxCreateDoubleMatrix(1_mwSize, 1_mwSize, mxREAL)
+  plhs(5) = mxCreateDoubleMatrix(n1, 1_mwSize, mxCOMPLEX)
+  plhs(6) = mxCreateDoubleMatrix(1_mwSize, 1_mwSize, mxREAL)
+
+  s => mxGetPr(plhs(1))
+  t => mxGetPr(plhs(2))
+  sdim => mxGetPr(plhs(4))
 #if MX_HAS_INTERLEAVED_COMPLEX
-  gev => mxGetComplexDoubles(plhs(6))
+  gev => mxGetComplexDoubles(plhs(5))
 #else
-  gev_r => mxGetPr(plhs(6))
-  gev_i => mxGetPi(plhs(6))
+  gev_r => mxGetPr(plhs(5))
+  gev_i => mxGetPi(plhs(5))
 #endif
-  info => mxGetPr(plhs(7))
-  z => mxGetPr(plhs(4))
+  info => mxGetPr(plhs(6))
+  z => mxGetPr(plhs(3))
 #if defined(MATLAB_MEX_FILE) && MATLAB_VERSION < 0x0904
   ! The left Schur vectors (VSL) are normally not computed, since JOBVSL="N".
   ! But old MKL versions (at least the one shipped with MATLAB R2009b/7.9,
@@ -162,6 +162,4 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   ! error number (only if no other error)
   if (any(abs(alpha_r) <= zhreshold .and. abs(beta) <= zhreshold) .and. info_bl == 0) &
        info = 30
-
-  plhs(1) = mxCreateDoubleScalar(0._real64)
 end subroutine mexFunction
diff --git a/mex/sources/ms-sbvar/mex_top_level.cc b/mex/sources/ms-sbvar/mex_top_level.cc
index e4060c2be946b8e0206bad969572a0a3a2b4614c..cee4610d1238a19faca767c20ba7b2ce42fbce0d 100644
--- a/mex/sources/ms-sbvar/mex_top_level.cc
+++ b/mex/sources/ms-sbvar/mex_top_level.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2010-2017 Dynare Team
+ * Copyright © 2010-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -42,8 +42,8 @@ mexFunction(int nlhs, mxArray *plhs[],
   /*
    * Check args
    */
-  if (nrhs != 1 || !mxIsChar(prhs[0]) || nlhs != 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("Error in MS-SBVAR MEX file: this function takes 1 string input argument and returns 1 output argument.");
+  if (nrhs != 1 || !mxIsChar(prhs[0]) || nlhs != 0)
+    mexErrMsgTxt("Error in MS-SBVAR MEX file: this function takes 1 string input argument and returns no output argument.");
 
   /*
    * Allocate memory
@@ -52,24 +52,24 @@ mexFunction(int nlhs, mxArray *plhs[],
   argument = static_cast<char *>(mxCalloc(mxGetN(prhs[0])+1, sizeof(char)));
   args = static_cast<char **>(mxCalloc(maxnargs, sizeof(char *)));
   if (argument == NULL || args == NULL)
-    DYN_MEX_FUNC_ERR_MSG_TXT("Error in MS-SBVAR MEX file: could not allocate memory. (1)");
+    mexErrMsgTxt("Error in MS-SBVAR MEX file: could not allocate memory. (1)");
 
   /*
    * Create argument string from prhs and parse to create args / nargs
    */
   if (!(args[nargs] = static_cast<char *>(mxCalloc(strlen(mainarg)+1, sizeof(char)))))
-    DYN_MEX_FUNC_ERR_MSG_TXT("Error in MS-SBVAR MEX file: could not allocate memory. (2)");
+    mexErrMsgTxt("Error in MS-SBVAR MEX file: could not allocate memory. (2)");
 
   strncpy(args[nargs++], mainarg, strlen(mainarg));
 
   if (mxGetString(prhs[0], argument, mxGetN(prhs[0])+1))
-    DYN_MEX_FUNC_ERR_MSG_TXT("Error in MS-SBVAR MEX file: error using mxGetString.\n");
+    mexErrMsgTxt("Error in MS-SBVAR MEX file: error using mxGetString.\n");
 
   beginarg = &argument[0];
   while ((n = strcspn(beginarg, " ")))
     {
       if (!(args[nargs] = static_cast<char *>(mxCalloc(n+1, sizeof(char)))))
-        DYN_MEX_FUNC_ERR_MSG_TXT("Error in MS-SBVAR MEX file: could not allocate memory. (3)");
+        mexErrMsgTxt("Error in MS-SBVAR MEX file: could not allocate memory. (3)");
 
       strncpy(args[nargs++], beginarg, n);
       beginarg += (isspace(beginarg[n]) || isblank(beginarg[n]) ? ++n : n);
@@ -85,7 +85,7 @@ mexFunction(int nlhs, mxArray *plhs[],
     }
   catch (const char *str)
     {
-      DYN_MEX_FUNC_ERR_MSG_TXT(str);
+      mexErrMsgTxt(str);
     }
 
   /*
@@ -94,6 +94,4 @@ mexFunction(int nlhs, mxArray *plhs[],
   for (n = 0; n < nargs; n++)
     mxFree(args[n]);
   mxFree(args);
-
-  plhs[0] = mxCreateDoubleScalar(0);
 }
diff --git a/mex/sources/sobol/qmc_sequence.cc b/mex/sources/sobol/qmc_sequence.cc
index 51b19b2e0fcef536b5edf6138dd33ec6b0603889..21a03dc1f9bf5cea6a4914f363a31f0158acadbf 100644
--- a/mex/sources/sobol/qmc_sequence.cc
+++ b/mex/sources/sobol/qmc_sequence.cc
@@ -1,7 +1,7 @@
 /*
 ** Computes Quasi Monte-Carlo sequence.
 **
-** Copyright © 2010-2017 Dynare Team
+** Copyright © 2010-2020 Dynare Team
 **
 ** This file is part of Dynare (can be used outside Dynare).
 **
@@ -52,23 +52,23 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   ** Check the number of input and output arguments.
   */
   if (nrhs < 3 || nrhs > 5)
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Five, four or three input arguments are required!");
+    mexErrMsgTxt("qmc_sequence:: Five, four or three input arguments are required!");
 
   if (nlhs == 0)
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: At least one output argument is required!");
+    mexErrMsgTxt("qmc_sequence:: At least one output argument is required!");
 
   /*
   ** Test the first input argument and assign it to dimension.
   */
   if (!mxIsNumeric(prhs[0]))
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: First input (dimension) has to be a positive integer!");
+    mexErrMsgTxt("qmc_sequence:: First input (dimension) has to be a positive integer!");
 
   int dimension = static_cast<int>(mxGetScalar(prhs[0]));
   /*
   ** Test the second input argument and assign it to seed.
   */
   if (!(mxIsNumeric(prhs[1]) && mxIsClass(prhs[1], "int64")))
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Second input (seed) has to be an integer [int64]!");
+    mexErrMsgTxt("qmc_sequence:: Second input (seed) has to be an integer [int64]!");
 
   int64_T seed = static_cast<int64_T>(mxGetScalar(prhs[1]));
   /*
@@ -83,22 +83,22 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     error_flag_3 = 1;
 
   if (error_flag_3 == 1)
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Third input (type of QMC sequence) has to be an integer equal to 0, 1 or 2!");
+    mexErrMsgTxt("qmc_sequence:: Third input (type of QMC sequence) has to be an integer equal to 0, 1 or 2!");
 
   /*
   ** Test dimension ≥ 2 when type==2
   */
   if (type == 2 && dimension < 2)
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: First input (dimension) has to be greater than 1 for a uniform QMC on an hypershere!");
+    mexErrMsgTxt("qmc_sequence:: First input (dimension) has to be greater than 1 for a uniform QMC on an hypershere!");
 
   else if (dimension > DIM_MAX)
-    DYN_MEX_FUNC_ERR_MSG_TXT(("qmc_sequence:: First input (dimension) has to be smaller than " + to_string(DIM_MAX) + " !").c_str());
+    mexErrMsgTxt(("qmc_sequence:: First input (dimension) has to be smaller than " + to_string(DIM_MAX) + " !").c_str());
 
   /*
   ** Test the optional fourth input argument and assign it to sequence_size.
   */
   if (nrhs > 3 && !mxIsNumeric(prhs[3]))
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: Fourth input (qmc sequence size) has to be a positive integer!");
+    mexErrMsgTxt("qmc_sequence:: Fourth input (qmc sequence size) has to be a positive integer!");
 
   int sequence_size;
   if (nrhs > 3)
@@ -110,17 +110,17 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   ** Test the optional fifth input argument and assign it to lower_and_upper_bounds.
   */
   if (nrhs > 4 && type == 0 && mxGetN(prhs[4]) != 2) // Sequence of uniformly distributed numbers in an hypercube
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be an array with two columns!");
+    mexErrMsgTxt("qmc_sequence:: The fifth input argument must be an array with two columns!");
 
   if (nrhs > 4 && type == 0 && static_cast<int>(mxGetM(prhs[4])) != dimension)
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fourth input argument must be an array with a number of lines equal to dimension (first input argument)!");
+    mexErrMsgTxt("qmc_sequence:: The fourth input argument must be an array with a number of lines equal to dimension (first input argument)!");
 
   if (nrhs > 4 && type == 1 && !(static_cast<int>(mxGetN(prhs[4])) == dimension
                                  && static_cast<int>(mxGetM(prhs[4])) == dimension)) // Sequence of normally distributed numbers
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be a squared matrix (whose dimension is given by the first input argument)!");
+    mexErrMsgTxt("qmc_sequence:: The fifth input argument must be a squared matrix (whose dimension is given by the first input argument)!");
 
   if (nrhs > 4 && type == 2 && !(mxGetN(prhs[4]) == 1 && mxGetM(prhs[4]) == 1)) // Sequence of uniformly distributed numbers on a hypershere
-    DYN_MEX_FUNC_ERR_MSG_TXT("qmc_sequence:: The fifth input argument must be a positive scalar!");
+    mexErrMsgTxt("qmc_sequence:: The fifth input argument must be a positive scalar!");
 
   const double *lower_bounds = nullptr, *upper_bounds = nullptr;
   int unit_hypercube_flag = 1;
diff --git a/tests/kalman_steady_state/test1.m b/tests/kalman_steady_state/test1.m
index d928d2b374526c58605dd4043e56c3abae966210..895aeb9f69602210e9c2a833c39dea3842a230b8 100644
--- a/tests/kalman_steady_state/test1.m
+++ b/tests/kalman_steady_state/test1.m
@@ -17,8 +17,7 @@ for i=1:100
     Z(1,1) = 1;
     Z(2,3) = 1;
     Z(3,6) = 1;
-    [err P] = kalman_steady_state(transpose(T),QQ,transpose(Z),H);
-    mexErrCheck('kalman_steady_state',err);
+    P = kalman_steady_state(transpose(T),QQ,transpose(Z),H);
 end
 
 % Without measurment errors.
@@ -35,6 +34,5 @@ for i=1:100
     Z(1,1) = 1;
     Z(2,3) = 1;
     Z(3,6) = 1;
-    [err P] = kalman_steady_state(transpose(T),QQ,transpose(Z),H);
-    mexErrCheck('kalman_steady_state',err);
+    P = kalman_steady_state(transpose(T),QQ,transpose(Z),H);
 end
diff --git a/tests/kronecker/test_kron.m b/tests/kronecker/test_kron.m
index a25b206689211b92a33e1cff8bf923c3c8e489e3..27e80877e04251e7e206b9ad633cc904726a39fb 100644
--- a/tests/kronecker/test_kron.m
+++ b/tests/kronecker/test_kron.m
@@ -1,5 +1,5 @@
 function info = test_kron(test,number_of_threads)
-% Copyright (C) 2007-2010 Dynare Team
+% Copyright (C) 2007-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -52,14 +52,12 @@ if test == 1
     disp('')
     disp('Computation of A*kron(B,B) with the mex file (v1):')
     tic 
-    [err, D1] = sparse_hessian_times_B_kronecker_C(A,B,number_of_threads);
-    mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
+    D1 = sparse_hessian_times_B_kronecker_C(A,B,number_of_threads);
     toc
     disp('')
     disp('Computation of A*kron(B,B) with the mex file (v2):')
     tic
-    [err, D2] = sparse_hessian_times_B_kronecker_C(A,B,B,number_of_threads);
-    mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
+    D2 = sparse_hessian_times_B_kronecker_C(A,B,B,number_of_threads);
     toc
     disp('');
     disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
@@ -120,14 +118,12 @@ if test==2
     disp(' ')
     disp('Computation of A*kron(B,B) with the mex file (v1):')
     tic
-    [err, D1] = sparse_hessian_times_B_kronecker_C(hessian,zx,number_of_threads);
-    mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
+    D1 = sparse_hessian_times_B_kronecker_C(hessian,zx,number_of_threads);
     toc
     disp(' ')
     disp('Computation of A*kron(B,B) with the mex file (v2):')
     tic
-    [err, D2] = sparse_hessian_times_B_kronecker_C(hessian,zx,zx,number_of_threads);
-    mexErrCheck('sparse_hessian_times_B_kronecker_C', err);    
+    D2 = sparse_hessian_times_B_kronecker_C(hessian,zx,zx,number_of_threads);
     toc
     disp(' ');
     disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
@@ -167,8 +163,7 @@ if test==3
     disp('Test with full format matrix -- 1(a)')
     D1 = A*kron(B,C);
     tic
-    [err, D2] = A_times_B_kronecker_C(A,B,C,number_of_threads);
-    mexErrCheck('A_times_B_kronecker_C', err);
+    D2 = A_times_B_kronecker_C(A,B,C,number_of_threads);
     toc
     disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
     if max(max(abs(D1-D2)))>1e-10
@@ -177,8 +172,7 @@ if test==3
     disp('Test with full format matrix -- 1(b)')
     D1 = A*kron(B,B);
     tic
-    [err, D2] = A_times_B_kronecker_C(A,B,number_of_threads);
-    mexErrCheck('A_times_B_kronecker_C', err);
+    D2 = A_times_B_kronecker_C(A,B,number_of_threads);
     toc
     disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
     if max(max(abs(D1-D2)))>1e-10
@@ -187,4 +181,4 @@ if test==3
     if ~(test3_1 && test3_2)
         info = 0;
     end
-end
\ No newline at end of file
+end