diff --git a/matlab/backward/backward_model_inversion.m b/matlab/backward/backward_model_inversion.m
index 8641b8d454b88c9d6ff5803804e0127845795c4b..0dd709e5b2823e402780dfec683e3e46621593ee 100644
--- a/matlab/backward/backward_model_inversion.m
+++ b/matlab/backward/backward_model_inversion.m
@@ -14,7 +14,7 @@ function [endogenousvariables, exogenousvariables] = backward_model_inversion(co
 %
 % REMARKS
 
-% Copyright (C) 2017 Dynare Team
+% Copyright © 2017-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -108,9 +108,9 @@ for t = 1:nobs(constraints)
     % values) and the free exogenous variables (initialized with 0).
     z = [Y(freeendogenousvariables_id,ity-1); zeros(nxfree, 1)];
     % Solves for z.
-    [z, failed] = dynare_solve(model_dtransf, z, DynareOptions, model_dynamic, ylag, ycur, X, DynareModel.params, DynareOutput.steady_state, itx, ModelInversion);
+    [z, failed, ~, ~, errorcode] = dynare_solve(model_dtransf, z, DynareOptions, model_dynamic, ylag, ycur, X, DynareModel.params, DynareOutput.steady_state, itx, ModelInversion);
     if failed
-        error('Enable to solve the system of equations!')
+        error('Nonlinear solver failed with errorcode=%i', errorcode)
     end
     % Update the matrix of exogenous variables.
     X(itx,freeinnovations_id) = z(nyfree+(1:nxfree));
diff --git a/matlab/backward/simul_backward_nonlinear_model_.m b/matlab/backward/simul_backward_nonlinear_model_.m
index bd87796bdd135fea6263e318c71767c69b2c3e2a..d84b605ed35957e97f105cd9ecaac5124a17bde6 100644
--- a/matlab/backward/simul_backward_nonlinear_model_.m
+++ b/matlab/backward/simul_backward_nonlinear_model_.m
@@ -56,17 +56,17 @@ for it = initialconditions.nobs+(1:samplesize)
     y = y_;                            % A good guess for the initial conditions is the previous values for the endogenous variables.
     try
         if ismember(DynareOptions.solve_algo, [12,14])
-            [DynareOutput.endo_simul(:,it), errorflag, ~, ~, exitflag] = ...
+            [DynareOutput.endo_simul(:,it), errorflag, ~, ~, errorcode] = ...
                 dynare_solve(model_dynamic_s, y, DynareOptions, ...
                              DynareModel.isloggedlhs, DynareModel.isauxdiffloggedrhs, DynareModel.endo_names, DynareModel.lhs, ...
                              model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
         else
-            [DynareOutput.endo_simul(:,it), errorflag, ~, ~, exitflag] = ...
+            [DynareOutput.endo_simul(:,it), errorflag, ~, ~, errorcode] = ...
                 dynare_solve(model_dynamic_s, y, DynareOptions, ...
                              model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
         end
         if errorflag
-            error()
+            error('Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
         end
     catch
         DynareOutput.endo_simul = DynareOutput.endo_simul(:, 1:it-1);
@@ -108,7 +108,7 @@ for it = initialconditions.nobs+(1:samplesize)
         %
         % Evaluate and check the residuals
         %
-        r = feval(model_dynamic, [ytm; ytm(iy1)], DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+        [r, J] = feval(model_dynamic, [ytm; ytm(iy1)], DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
         residuals_evaluating_to_nan = isnan(r);
         residuals_evaluating_to_inf = isinf(r);
         residuals_evaluating_to_complex = ~isreal(r);
@@ -130,34 +130,6 @@ for it = initialconditions.nobs+(1:samplesize)
             display_names_of_problematic_equations(DynareModel, residuals_evaluating_to_nan);
             skipline()
         end
-        %
-        % Display value of exitflag if available (fsolve, solve_algo=0, only)
-        %
-        if ~isnan(exitflag)
-            skipline()
-            switch exitflag
-              case 0
-                disp('Returned value for exitflag is 0 (maximum number of iterations or evaluation reached).')
-              case -1
-                if options.solve_algo==0
-                    disp('Returned value for exitflag is -1 (objective function stopped algorithm).')
-                elseif options.solve_algo==1
-                    disp('Returned value for exitflag is -1 (objective function cannot be evaluated at the initial guess).')
-                end
-              case -3
-                if options.solve_algo==0
-                    disp('Returned value for exitflag is -3 (Trust region radius became too small, trust-region-dogleg algorithm).')
-                elseif options.solve_algo==1
-                    disp('Returned value for exitflag is -3 (Convergence too slow).')
-                end
-              case -2
-                if options.solve_algo==0
-                    disp('Returned value for exitflag is -2.')
-                elseif options.solve_algo==1
-                    disp('Returned value for exitflag is -2 (Spurious convergence).')
-                end
-            end
-        end
         break
         % TODO Implement same checks with the jacobian matrix.
         % TODO Modify other solvers to return an exitflag.
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index ff2be72f97b6bc3fb15e109b8d6238b981a8ec93..985ddd47910ca22014165d352f75c6ab52fb05ad 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -1,4 +1,4 @@
-function [x, errorflag, fvec, fjac, exitflag] = dynare_solve(f, x, options, varargin)
+function [x, errorflag, fvec, fjac, errorcode] = dynare_solve(f, x, options, varargin)
 
 % Solves a nonlinear system of equations, f(x) = 0 with n unknowns
 % and n equations.
@@ -14,7 +14,13 @@ function [x, errorflag, fvec, fjac, exitflag] = dynare_solve(f, x, options, vara
 % - errorflag    [logical]          scalar, true iff the model can not be solved.
 % - fvec         [double]           n×1 vector, function value at x (f(x), used for debugging when errorflag is true).
 % - fjac         [double]           n×n matrix, Jacobian value at x (J(x), used for debugging when errorflag is true).
-% - exitflag     [integer]          scalar,
+% - errorcode    [integer]          scalar.
+%
+% REMARKS
+% Interpretation of the error code depends on the algorithm, except if value of errorcode is
+%
+%        -10  -> System of equation ill-behaved at the initial guess (Inf, Nans or complex numbers).
+%        -11  -> Initial guess is a solution of the system of equations.
 
 % Copyright © 2001-2022 Dynare Team
 %
@@ -33,8 +39,6 @@ function [x, errorflag, fvec, fjac, exitflag] = dynare_solve(f, x, options, vara
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-exitflag = nan;
-
 jacobian_flag = options.jacobian_flag; % true iff Jacobian is returned by f routine (as a second output argument).
 
 % Set tolerance parameter depending the the caller function.
@@ -49,7 +53,7 @@ else
         caller_file_name=stack(1).file;
     end
 end
-if strcmp(caller_file_name, 'solve_stacked_problem.m') || strcmp(caller_file_name, 'sim1_purely_backward.m') 
+if strcmp(caller_file_name, 'solve_stacked_problem.m') || strcmp(caller_file_name, 'sim1_purely_backward.m')
     tolf = options.dynatol.f;
     tolx = options.dynatol.x;
 else
@@ -63,7 +67,7 @@ else
     maxit = options.steady.maxit;
 end
 
-errorflag = false;
+errorflag = false; % Let's be optimistic!
 nn = size(x,1);
 
 % Keep a copy of the initial guess.
@@ -98,7 +102,7 @@ if jacobian_flag
     if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)))
         if ~ismember(options.solve_algo,[10,11]) && max(abs(fvec))< tolf
             % return if initial value solves the problem except if a mixed complementarity problem is to be solved (complementarity conditions may not be satisfied)
-            exitflag = -1;
+            errorcode = -11;
             return;
         end
         disp_verbose('Randomize initial guess...', options.verbosity)
@@ -134,7 +138,7 @@ else
     fjac = zeros(nn, nn);
     if ~ismember(options.solve_algo,[10,11]) && max(abs(fvec)) < tolf
         % return if initial value solves the problem except if a mixed complementarity problem is to be solved (complementarity conditions may not be satisfied)
-        exitflag = -1;
+        errorcode = -11;
         return;
     end
     wrong_initial_guess_flag = false;
@@ -170,6 +174,7 @@ end
 
 % Exit with error if no initial guess has been found.
 if wrong_initial_guess_flag
+    errorcode = -10;
     errorflag = true;
     x = x0;
     return
@@ -196,7 +201,7 @@ if options.solve_algo == 0
         options4fsolve.Jacobian = 'off';
     end
     if ~isoctave
-        [x, ~, exitflag] = fsolve(f, x, options4fsolve, arguments{:});
+        [x, ~, errorcode] = fsolve(f, x, options4fsolve, arguments{:});
     else
         % Under Octave, use a wrapper, since fsolve() does not have a 4th arg
         if ischar(f)
@@ -205,19 +210,20 @@ if options.solve_algo == 0
             f2 = f;
         end
         f = @(x) f2(x, arguments{:});
-        [x, ~, exitflag] = fsolve(f, x, options4fsolve);
+        [x, ~, errorcode] = fsolve(f, x, options4fsolve);
     end
-    if exitflag == 1
+    if errorcode==1
         errorflag = false;
-    elseif exitflag > 1
-        if ischar(f)
-            f2 = str2func(f);
-        else
-            f2 = f;
+    elseif errorcode>1
+        if ~isoctave
+            if ischar(f)
+                f2 = str2func(f);
+            else
+                f2 = f;
+            end
+            f = @(x) f2(x, arguments{:});
         end
-        f = @(x) f2(x, arguments{:});
-        fvec = feval(f, x);
-        if max(abs(fvec)) >= tolf
+        if max(abs(fvec)) > tolf
             errorflag = true;
         else
             errorflag = false;
@@ -225,10 +231,13 @@ if options.solve_algo == 0
     else
         errorflag = true;
     end
+    [fvec, fjac] = feval(f, x, arguments{:});
 elseif options.solve_algo==1
-    [x, errorflag, exitflag] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, arguments{:});
+    [x, errorflag, errorcode] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, arguments{:});
+    [fvec, fjac] = feval(f, x, arguments{:});
 elseif options.solve_algo==9
-    [x, errorflag, exitflag] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, arguments{:});
+    [x, errorflag, errorcode] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, arguments{:});
+    [fvec, fjac] = feval(f, x, arguments{:});
 elseif ismember(options.solve_algo, [2, 12, 4])
     if ismember(options.solve_algo, [2, 12])
         solver = @solve1;
@@ -304,11 +313,11 @@ elseif ismember(options.solve_algo, [2, 12, 4])
                 dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u with trust_region routine.', i);
             end
         end
-        [x, errorflag, exitflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
-                                          options.gstep, ...
-                                          tolf, options.solve_tolx, maxit, ...
-                                          options.trust_region_initial_step_bound_factor, ...
-                                          options.debug, arguments{:});
+        [x, errorflag, errorcode] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
+                                           options.gstep, ...
+                                           tolf, options.solve_tolx, maxit, ...
+                                           options.trust_region_initial_step_bound_factor, ...
+                                           options.debug, arguments{:});
         fre = true;
         if errorflag
             return
@@ -317,27 +326,34 @@ elseif ismember(options.solve_algo, [2, 12, 4])
     fvec = feval(f, x, arguments{:});
     if max(abs(fvec))>tolf
         disp_verbose('Call solver on the full nonlinear problem.',options.verbosity)
-        [x, errorflag, exitflag] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
-                                          options.gstep, tolf, options.solve_tolx, maxit, ...
-                                          options.trust_region_initial_step_bound_factor, ...
-                                          options.debug, arguments{:});
+        [x, errorflag, errorcode] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
+                                           options.gstep, tolf, options.solve_tolx, maxit, ...
+                                           options.trust_region_initial_step_bound_factor, ...
+                                           options.debug, arguments{:});
     end
+    [fvec, fjac] = feval(f, x, arguments{:});
 elseif options.solve_algo==3
     if jacobian_flag
-        [x, errorflag] = csolve(f, x, f, tolf, maxit, arguments{:});
+        [x, errorcode] = csolve(f, x, f, tolf, maxit, arguments{:});
     else
-        [x, errorflag] = csolve(f, x, [], tolf, maxit, arguments{:});
+        [x, errorcode] = csolve(f, x, [], tolf, maxit, arguments{:});
+    end
+    if errorcode==0
+        errorflag = false;
+    else
+        errorflag = true;
     end
     [fvec, fjac] = feval(f, x, arguments{:});
 elseif options.solve_algo==10
     % LMMCP
     olmmcp = options.lmmcp;
-    [x, ~, exitflag] = lmmcp(f, x, olmmcp.lb, olmmcp.ub, olmmcp, arguments{:});
-    if exitflag==1
+    [x, fvec, errorcode, ~, fjac] = lmmcp(f, x, olmmcp.lb, olmmcp.ub, olmmcp, arguments{:});
+    if errorcode==1
         errorflag = false;
     else
         errorflag = true;
     end
+    [fvec, fjac] = feval(f, x, arguments{:});
 elseif options.solve_algo == 11
     % PATH mixed complementary problem
     % PATH linear mixed complementary problem
@@ -355,6 +371,8 @@ elseif options.solve_algo == 11
     catch
         errorflag = true;
     end
+    errorcode = nan; % There is no error code for this algorithm, as PATH is closed source it is unlikely we can fix that.
+    [fvec, fjac] = feval(f, x, arguments{:});
 elseif ismember(options.solve_algo, [13, 14])
     if ~jacobian_flag
         error('DYNARE_SOLVE: option solve_algo=13|14 needs computed Jacobian')
@@ -366,7 +384,7 @@ elseif ismember(options.solve_algo, [13, 14])
         auxstruct.isloggedlhs = isloggedlhs;
         auxstruct.isauxdiffloggedrhs = isauxdiffloggedrhs;
     end
-    [x, errorflag, exitflag] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, auxstruct, arguments{:});
+    [x, errorflag, errorcode] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, auxstruct, arguments{:});
     [fvec, fjac] = feval(f, x, arguments{:});
 else
     error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
index 3a9c0d88f05e7e5f11674cdeb1b673bc40f0f7fd..d3f7ec78d275c53aaa25b34a803d2b005946fec3 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
@@ -1,6 +1,6 @@
-function [flag,endo_simul,err,y] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,Options,pfm,order,varargin)
+function [errorflag, endo_simul, errorcode, y] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, Options, pfm, order, varargin)
 
-% Copyright (C) 2012-2017 Dynare Team
+% Copyright © 2012-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -22,9 +22,9 @@ if nargin < 6
 else
     homotopy_parameter = varargin{1};
 end
+
 flag = 0;
 err = 0;
-stop = 0;
 
 EpOptions = Options.ep;
 
@@ -143,7 +143,7 @@ Options.steady.maxit = 100;
 y = repmat(steady_state,block_nbr,1);
 Options.solve_algo = Options.ep.solve_algo;
 Options.steady.maxit = Options.ep.maxit;
-[y,info] = dynare_solve(@ep_problem_2,y,Options,exo_simul,pfm);
+[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_2,y,Options,exo_simul,pfm);
 if info
     flag = 1;
     err = info;
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index d369d9586bcffa13deb41acd287b4098fe382eb4..651363dc5cebd4e338a8b4db16554591e769e8b9 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-2021 Dynare Team
+% Copyright © 2001-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -262,17 +262,19 @@ elseif ~options.bytecode && ~options.block
     if ~options.linear
         % non linear model
         static_model = str2func([M.fname '.static']);
-        [ys,check] = dynare_solve(@static_problem,...
-                                  ys_init,...
-                                  options, exo_ss, params,...
-                                  M.endo_nbr,...
-                                  static_model);
+        [ys, check] = dynare_solve(@static_problem, ...
+                                   ys_init, ...
+                                   options, exo_ss, params, ...
+                                   M.endo_nbr, ...
+                                   static_model);
         if check && options.debug
-            [ys,check,fvec,fjac] = dynare_solve(@static_problem,...
-                                                ys_init,...
-                                                options, exo_ss, params,...
-                                                M.endo_nbr,...
-                                                static_model);
+            [ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ...
+                                                              ys_init, ...
+                                                              options, exo_ss, params, ...
+                                                              M.endo_nbr, ...
+                                                              static_model);
+            dprintf('Nonlinear solver routine returned errorcode=%i.', errorcode)
+            skipline()
             [infrow,infcol]=find(isinf(fjac) | isnan(fjac));
             if ~isempty(infrow)
                 fprintf('\nSTEADY:  The Jacobian at the initial values contains Inf or NaN. The problem arises from: \n')
diff --git a/matlab/evaluate_steady_state_file.m b/matlab/evaluate_steady_state_file.m
index f2c9ec6b42fe6cc38a8ad7d4364c7704379e9518..4a804962f7f54dc56ba562642018f02e9bfe6792 100644
--- a/matlab/evaluate_steady_state_file.m
+++ b/matlab/evaluate_steady_state_file.m
@@ -19,7 +19,7 @@ function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options,
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2001-2019 Dynare Team
+% Copyright © 2001-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -127,5 +127,5 @@ elseif ~isempty(options.steadystate_partial)
     for i = 1:nov
         indv(i) = strmatch(ssvar(i), M.endo_names, 'exact');
     end
-    [ys,~] = dynare_solve('restricted_steadystate', ys(indv), options, exo_ss,indv);
+    ys = dynare_solve('restricted_steadystate', ys(indv), options, exo_ss,indv);
 end
diff --git a/matlab/perfect-foresight-models/sim1_purely_backward.m b/matlab/perfect-foresight-models/sim1_purely_backward.m
index 0641b55a12333c72f7f86db3c3dade08065f3b68..81c5e9edc342170836130c4cebe179552bbbee4b 100644
--- a/matlab/perfect-foresight-models/sim1_purely_backward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_backward.m
@@ -2,7 +2,7 @@ function [endogenousvariables, info] = sim1_purely_backward(endogenousvariables,
 
 % Performs deterministic simulation of a purely backward model
 
-% Copyright © 2012-2021 Dynare Team
+% Copyright © 2012-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -39,14 +39,17 @@ for it = M.maximum_lag + (1:options.periods)
     y = endogenousvariables(:,it-1);        % Values at previous period, also used as guess value for current period
     ylag = y(iyb);
     if ismember(options.solve_algo, [12,14])
-        [tmp, check] = dynare_solve(dynamicmodel_s, y, options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
-                                    dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
+        [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
+                                                     dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
     else
-        [tmp, check] = dynare_solve(dynamicmodel_s, y, options, ...
-                                    dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
+        [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, options, ...
+                                                     dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
     end
     if check
         info.status = false;
+        if options.debug
+            dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i. in period %i', errorcode, it)
+        end
     end
     endogenousvariables(:,it) = tmp;
 end
\ No newline at end of file
diff --git a/matlab/perfect-foresight-models/sim1_purely_static.m b/matlab/perfect-foresight-models/sim1_purely_static.m
index b1bf5591f112876a7e4f3f381e88c10fc9efc040..373896cff9e54db8a05e2f7c03562e35b6684767 100644
--- a/matlab/perfect-foresight-models/sim1_purely_static.m
+++ b/matlab/perfect-foresight-models/sim1_purely_static.m
@@ -2,7 +2,7 @@ function [endogenousvariables, info] = sim1_purely_static(endogenousvariables, e
 
 % Performs deterministic simulation of a purely static model
 
-% Copyright © 2021 Dynare Team
+% Copyright © 2021-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,14 +36,17 @@ y = endogenousvariables(:,1);
 
 for it = 1:options.periods
     if ismember(options.solve_algo, [12,14])
-        [tmp, check] = dynare_solve(dynamicmodel_s, y, options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
-                                    dynamicmodel, exogenousvariables, M.params, steadystate, it);
+        [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
+                                                     dynamicmodel, exogenousvariables, M.params, steadystate, it);
     else
-        [tmp, check] = dynare_solve(dynamicmodel_s, y, options, ...
-                                    dynamicmodel, exogenousvariables, M.params, steadystate, it);
+        [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, options, ...
+                                                     dynamicmodel, exogenousvariables, M.params, steadystate, it);
     end
     if check
         info.status = false;
+        if options.debug
+            dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+        end
     end
     endogenousvariables(:,it) = tmp;
     y = endogenousvariables(:,it);
diff --git a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
index ca54b95310488f41cd0833057f92572212aa9651..2a30699de35bbf56fc5d29bc7a31d09e3c57855b 100644
--- a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
@@ -1,6 +1,6 @@
 function [endogenousvariables, info] = solve_stacked_linear_problem(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
 
-% Copyright (C) 2015-2019 Dynare Team
+% Copyright © 2015-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -41,12 +41,12 @@ jendo = transpose(1:nd);
 z = bsxfun(@minus, z, steadystate_y);
 x = bsxfun(@minus, exogenousvariables, steadystate_x');
 
-[y, check] = dynare_solve(@linear_perfect_foresight_problem,z(:), options, ...
-                          jacobian, y0-steadystate_y, yT-steadystate_y, ...
-                          x, M.params, steadystate_y, ...
-                          M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
-                          i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
-                          jendo, jexog);
+[y, check, ~, ~, errorcode] = dynare_solve(@linear_perfect_foresight_problem,z(:), options, ...
+                                           jacobian, y0-steadystate_y, yT-steadystate_y, ...
+                                           x, M.params, steadystate_y, ...
+                                           M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
+                                           i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
+                                           jendo, jexog);
 
 if all(imag(y)<.1*options.dynatol.x)
     if ~isreal(y)
@@ -60,6 +60,9 @@ endogenousvariables = [y0 bsxfun(@plus,reshape(y,M.endo_nbr,options.periods), st
 
 if check
     info.status = false;
+    if options.debug
+        dprintf('solve_stacked_linear_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
+    end
 else
     info.status = true;
 end
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
index 53054c4f5f94e9040302ab7aab3a098671d86c76..484d7f9cf6d580321185c33a7627b0b02a07c7b8 100644
--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -13,7 +13,7 @@ function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables
 % - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
 % - info                [struct] contains informations about the results.
 
-% Copyright (C) 2015-2019 Dynare Team
+% Copyright © 2015-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -46,14 +46,14 @@ if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity
         options.mcppath.lb = repmat(lb,options.periods,1);
         options.mcppath.ub = repmat(ub,options.periods,1);
     end
-    [y, check] = dynare_solve(@perfect_foresight_mcp_problem,z(:),options, ...
-                              dynamicmodel, y0, yT, ...
-                              exogenousvariables, M.params, steadystate, ...
-                              M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
-                              i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
-                              eq_index);
+    [y, check, ~, ~, errorcode] = dynare_solve(@perfect_foresight_mcp_problem,z(:),options, ...
+                                               dynamicmodel, y0, yT, ...
+                                               exogenousvariables, M.params, steadystate, ...
+                                               M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
+                                               i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
+                                               eq_index);
 else
-    [y, check] = dynare_solve(@perfect_foresight_problem,z(:), options, y0, yT, exogenousvariables, M.params, steadystate, options.periods, M, options);
+    [y, check, ~, ~, errorcode] = dynare_solve(@perfect_foresight_problem,z(:), options, y0, yT, exogenousvariables, M.params, steadystate, options.periods, M, options);
 end
 
 if all(imag(y)<.1*options.dynatol.x)
@@ -68,6 +68,9 @@ endogenousvariables(:, M.maximum_lag+(1:options.periods)) = reshape(y, M.endo_nb
 
 if check
     info.status = false;
+    if options.debug
+        dprintf('solve_stacked_problem: Nonlinear solver routine failed with errorcode=%i.', errorcode)
+    end
 else
     info.status = true;
 end
diff --git a/matlab/simul_static_model.m b/matlab/simul_static_model.m
index 97c3726b10f8632eaebd4bcf66ba508bff68dabb..44e8d09e0a101fe548db03012d55045c4425491b 100644
--- a/matlab/simul_static_model.m
+++ b/matlab/simul_static_model.m
@@ -15,7 +15,7 @@ function simulation = simul_static_model(samplesize, innovations)
 % [2] The last input argument is not mandatory. If absent we use random draws and rescale them with the informations provided
 %     through the shocks block.
 
-% Copyright (C) 2019 Dynare Team
+% Copyright © 2019-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -88,9 +88,9 @@ staticmodel = str2fun(sprintf('%s.static', M_.fname));
 % Simulations (call a Newton-like algorithm for each period).
 for t=1:samplesize
     y = zeros(M_.endo_nbr, 1);
-    [oo_.endo_simul(:,t), info] = dynare_solve(staticmodel, y, options_, oo_.exo_simul(t,:), M_.params);
+    [oo_.endo_simul(:,t), info, ~, ~, errorcode] = dynare_solve(staticmodel, y, options_, oo_.exo_simul(t,:), M_.params);
     if info
-        error('Newton failed!')
+        error('Newton failed (with error code %i).', errorcode)
     end
 end
 
diff --git a/matlab/solve1.m b/matlab/solve1.m
index 7637f057f03951e6b2b27159f350eebd111e7d8a..fa6a3262fdb72f7e77230ceee43f184117c3f6df 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -1,4 +1,4 @@
-function [x, errorflag, exitflag] = solve1(func, x, j1, j2, jacobian_flag, gstep, tolf, tolx, maxit, fake, debug, varargin)
+function [x, errorflag, errorcode] = solve1(func, x, j1, j2, jacobian_flag, gstep, tolf, tolx, maxit, fake, debug, varargin)
 
 % Solves systems of non linear equations of several variables
 %
@@ -61,7 +61,7 @@ idCpx = ~isreal(fvec);
 if any(idInf)
     disp('SOLVE1: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number:')
     disp(j1(idInf)')
-    exitflag = -1;
+    errorcode = 0;
     errorflag = true;
     return
 end
@@ -69,7 +69,7 @@ end
 if any(idNan)
     disp('SOLVE1: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a nan:')
     disp(j1(idNan)')
-    exitflag = -1;
+    errorcode = 0;
     errorflag = true;
     return
 end
@@ -77,7 +77,7 @@ end
 if any(idNan)
     disp('SOLVE1: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a complex number:')
     disp(j1(idCpx)')
-    exitflag = -1;
+    errorcode = 0;
     errorflag = true;
     return
 end
@@ -86,7 +86,7 @@ f = 0.5*(fvec'*fvec);
 
 if max(abs(fvec))<tolf*tolf
     % Initial guess is a solution
-    exitflag = 1;
+    errorcode = -1;
     return
 end
 
@@ -140,7 +140,7 @@ for its = 1:maxit
         den = max([f;0.5*nn]) ;
         if max(abs(g).*max([abs(x(j2)') ones(1,nn)])')/den < tolmin
             if max(abs(x(j2)-xold(j2))./max([abs(x(j2)') ones(1,nn)])') < tolx
-                exitflag = -3;
+                errorcode = 3;
                 if nargout<3
                     skipline()
                     dprintf('SOLVE: Iteration %s', num2str(its))
@@ -150,7 +150,7 @@ for its = 1:maxit
                 return
             end
         else
-            exitflag = -2;
+            errorcode = 4;
             if nargout<3
                 skipline()
                 dprintf('SOLVE: Iteration %s', num2str(its))
@@ -160,18 +160,19 @@ for its = 1:maxit
             return
         end
     elseif max(abs(fvec)) < tolf
-        exitflag = 1;
+        errorcode = 1;
         return
     end
 end
 
 errorflag = true;
-exitflag = 0;
+errorcode = 2;
 
 if nargout<3
     skipline()
     disp('SOLVE: maxit has been reached')
 end
+
 % 01/14/01 MJ lnsearch is now a separate function
 % 01/16/01 MJ added varargin to function evaluation
 % 04/13/01 MJ added test  f < tolf !!
diff --git a/matlab/static_model_inversion.m b/matlab/static_model_inversion.m
index 9993d0503cbba2a1cba0846f9a6374c239c3223e..c2b8d80e2f4990c88f69dc40a64d1f99b15a814d 100644
--- a/matlab/static_model_inversion.m
+++ b/matlab/static_model_inversion.m
@@ -13,7 +13,7 @@ function [endogenousvariables, exogenousvariables] = static_model_inversion(cons
 %
 % REMARKS
 
-% Copyright (C) 2019 Dynare Team
+% Copyright © 2019-2022 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -103,9 +103,9 @@ for t = 1:nobs(constraints)
     % values) and the free exogenous variables (initialized with 0).
     z = [Y(freeendogenousvariables_id,ity); zeros(nxfree, 1)];
     % Solves for z.
-    [z, failed] = dynare_solve(model_stransf, z, DynareOptions, model_dynamic, ycur, X(itx, :), DynareModel.params, ModelInversion);
-    if failed
-        error('Enable to solve the system of equations!')
+    [z, errorflag, ~, ~, errorcode] = dynare_solve(model_stransf, z, DynareOptions, model_dynamic, ycur, X(itx, :), DynareModel.params, ModelInversion);
+    if errorflag
+        error('Enable to solve the system of equations (with error code %i).', errorcode)
     end
     % Update the matrix of exogenous variables.
     X(itx,freeinnovations_id) = z(nyfree+(1:nxfree));