From 7e770f69e791e705ab535126e0157130dc9d40c6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 10 Jan 2020 17:55:57 +0100
Subject: [PATCH] Remove workaround for errors in MEX files

Because at some point throwing exceptions from MEX files (with mexErrMsgTxt())
was not working under Windows 64-bit, we had designed a workaround to avoid
using exceptions.

Most MEX files were returning an error code as their first (or sometimes last)
argument, and that code would have to be checked from the MATLAB code.

Since this workaround is no longer needed, this commit removes it. As a
consequence, the interface of many MEX files is modified.

For some background, see https://www.dynare.org/pipermail/dev/2010-September/000895.html
---
 matlab/DsgeSmoother.m                         |   5 +-
 matlab/block_bytecode_mfs_steadystate.m       |   4 +-
 matlab/bytecode_steadystate.m                 |   4 +-
 matlab/dr_block.m                             |  10 +-
 matlab/dsge_likelihood.m                      |  27 +-
 matlab/dyn_first_order_solver.m               |   5 +-
 matlab/dyn_ramsey_static.m                    |  10 +-
 matlab/dyn_risky_steadystate_solver.m         |   3 +-
 matlab/dyn_second_order_solver.m              |  23 +-
 matlab/dynare_solve_block_or_bytecode.m       |  19 +-
 matlab/ep/extended_path_core.m                |  12 +-
 matlab/evaluate_planner_objective.m           |  32 +-
 matlab/evaluate_static_model.m                |   7 +-
 matlab/evaluate_steady_state.m                |   5 +-
 matlab/get_perturbation_params_derivs.m       |  44 +--
 matlab/k_order_pert.m                         |   7 +-
 matlab/mex/k_order_perturbation.m             |   5 +-
 matlab/mexErrCheck.m                          |  42 ---
 matlab/missing/mex/gensylv/gensylv.m          |   7 +-
 .../mex/kronecker/A_times_B_kronecker_C.m     |   6 +-
 .../sparse_hessian_times_B_kronecker_C.m      |  13 +-
 matlab/missing/mex/mjdgges/mjdgges.m          |   7 +-
 matlab/model_diagnostics.m                    |  18 +-
 matlab/ms-sbvar/ms_compute_mdd.m              |   5 +-
 matlab/ms-sbvar/ms_compute_probabilities.m    |   5 +-
 matlab/ms-sbvar/ms_estimation.m               |   5 +-
 matlab/ms-sbvar/ms_forecast.m                 |   5 +-
 matlab/ms-sbvar/ms_irf.m                      |   5 +-
 matlab/ms-sbvar/ms_sbvar_setup.m              |   5 +-
 matlab/ms-sbvar/ms_simulation.m               |   5 +-
 matlab/ms-sbvar/ms_variance_decomposition.m   |   5 +-
 .../add_auxiliary_variables_to_steadystate.m  |   8 +-
 .../det_cond_forecast.m                       |  19 +-
 .../perfect_foresight_solver_core.m           |  40 +--
 matlab/resid.m                                |   5 +-
 matlab/simult_.m                              |  54 +--
 matlab/solve_perfect_foresight_model.m        |  12 +-
 matlab/stochastic_solvers.m                   |  10 +-
 .../block_kalman_filter.cc                    |  30 +-
 mex/sources/bytecode/bytecode.cc              | 333 ++++++++----------
 mex/sources/bytecode/testing/mex_interface.hh |   9 +-
 mex/sources/dynare_simul_/dynare_simul_.cc    |  37 +-
 mex/sources/dynmex.h                          |  16 +-
 mex/sources/gensylv/gensylv.cc                |  29 +-
 .../k_order_perturbation.cc                   |  85 +++--
 .../kalman_steady_state.cc                    |  47 ++-
 .../kronecker/A_times_B_kronecker_C.cc        |  15 +-
 .../sparse_hessian_times_B_kronecker_C.cc     |  17 +-
 mex/sources/mjdgges/mjdgges.F08               |  32 +-
 mex/sources/ms-sbvar/mex_top_level.cc         |  18 +-
 mex/sources/sobol/qmc_sequence.cc             |  26 +-
 tests/kalman_steady_state/test1.m             |   6 +-
 tests/kronecker/test_kron.m                   |  22 +-
 53 files changed, 527 insertions(+), 698 deletions(-)
 delete mode 100644 matlab/mexErrCheck.m

diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 793afa5616..ee1546af50 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 a34b15c197..ba0bef8617 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 4c3e2f46f9..dff35373f4 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 0328c1048c..d65248b774 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 40a7203c43..ddcdf33f40 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 63cb4f9da5..101a33700c 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 e94deb5044..6ea85e7d6e 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 d8fe2ede16..93c75da4dc 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 4e3151497c..bd5c37a9fa 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 c995c93c7d..f991f54f53 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 1768c0bbd6..eb8b2454dd 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 198d65ccda..63ec0ab6a7 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 d15005aa20..fbcfa235cd 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 12507a2b46..8fc5826c87 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 64c5a2e346..52b5f9cecf 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 7627970cec..beaa621b68 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 2f4b1eec96..15e34a7ce8 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 708d4e5057..0000000000
--- 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 2920aca4d2..462e7b1ea0 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 7a8d09498e..db615145a1 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 fd6059f840..7857e2571e 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 6bd345a570..11d51d62f7 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 5696f77269..70b39fe707 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 499baf6a52..43372a6e0a 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 07b951f65a..9a5b91de89 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 495a483eed..2a8510a2a4 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 e8949c9ae2..b6e25aa847 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 47d1a8490b..d06a1149c0 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 a195cc4913..062eaca728 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 f7832e035f..d311733b80 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 1fe73ddbc4..746facf2ca 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 e1ffe4f196..35bd84bb33 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 98875fb917..cdcbe24d77 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 85b420cd4b..267545d692 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 716c9c5b36..1f3743a072 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 e3ec1ec317..508af26e43 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 cedf3e0604..93c5fb26fb 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 5aa9920b41..87d88bacf1 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 0f484b04ff..0ad1f60c80 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 ee81f0843c..31efcd311c 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 9e89891893..b8b0dbfb0a 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 58a6396092..faaf146b03 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 52a71bef5a..4bb0ca7391 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 15a980ce7b..c71f295f92 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 0a44e82062..4ae1151db7 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 0a0db31818..14b61015b3 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 071122e3f9..95d0f72680 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 93265998cc..5eaa7eee59 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 dbc5da9d12..5209c01d1d 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 e4060c2be9..cee4610d12 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 51b19b2e0f..21a03dc1f9 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 d928d2b374..895aeb9f69 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 a25b206689..27e80877e0 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
-- 
GitLab