diff --git a/matlab/backward/backward_model_forecast.m b/matlab/backward/backward_model_forecast.m
index 87c6c424f0682c4d837df9fdf172d6547542f3f0..8f4ec8c0a5a6782a64761a04d46c26fa39fb40c4 100644
--- a/matlab/backward/backward_model_forecast.m
+++ b/matlab/backward/backward_model_forecast.m
@@ -10,7 +10,7 @@ function forecasts = backward_model_forecast(initialcondition, listofvariables,
 % OUTPUTS
 % - forecast            [dseries]
 
-% Copyright © 2017-2018 Dynare Team
+% Copyright © 2017-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -91,10 +91,15 @@ end
 
 % Compute forecast without shock
 if options_.linear
-    ysim__0 = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic);
+    [ysim__0, errorflag] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic);
 else
-    ysim__0 = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic);
+    [ysim__0, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic);
 end
+
+if errorflag
+    error('Simulation failed.')
+end
+
 forecasts.pointforecast = dseries(transpose(ysim__0(idy,:)), initialcondition.init, listofvariables);
 
 % Set first period of forecast
@@ -106,9 +111,12 @@ if withuncertainty
     for i=1:B
         innovations = transpose(sigma*randn(M_.exo_nbr, periods));
         if options_.linear
-            [ysim__, xsim__] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic);
+            [ysim__, xsim__, errorflag] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic);
         else
-            [ysim__, xsim__] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic);
+            [ysim__, xsim__, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic);
+        end
+        if errorflag
+            error('Simulation failed.')
         end
         ArrayOfForecasts(:,:,i) = ysim__(idy,:);
     end
@@ -120,4 +128,4 @@ if withuncertainty
     ArrayOfForecasts = sort(ArrayOfForecasts, 3);
     forecasts.lb = dseries(transpose(ArrayOfForecasts(:,:,round(0.025*B))), initialcondition.init, listofvariables);
     forecasts.ub = dseries(transpose(ArrayOfForecasts(:,:,round(0.975*B))), initialcondition.init, listofvariables);
-end
\ No newline at end of file
+end
diff --git a/matlab/backward/backward_model_irf.m b/matlab/backward/backward_model_irf.m
index 93b4fae717d1eb44fbec8bc9f328ec60df180786..ab1418e7802a96fe90202bfe6d537eea8ba019e6 100644
--- a/matlab/backward/backward_model_irf.m
+++ b/matlab/backward/backward_model_irf.m
@@ -19,7 +19,7 @@ function [deviations, baseline, irfs] = backward_model_irf(initialcondition, inn
 %   argument.
 % - If second argument is not empty, periods must not be greater than innovationbaseline.nobs.
 
-% Copyright © 2017-2020 Dynare Team
+% Copyright © 2017-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -163,9 +163,13 @@ irfs = struct();
 % Baseline paths (get transition paths induced by the initial condition and
 % baseline innovations).
 if options_.linear
-    ysim__0 = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, Innovations, nx, ny1, iy1, jdx, model_dynamic);
+    [ysim__0, errorflag] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, Innovations, nx, ny1, iy1, jdx, model_dynamic);
 else
-    ysim__0 = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, Innovations, iy1, model_dynamic);
+    [ysim__0, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, Innovations, iy1, model_dynamic);
+end
+
+if errorflag
+    error('Simulation failed. Cannot simulate baseline.')
 end
 
 % Transform the endogenous variables.
@@ -198,9 +202,13 @@ for i=1:length(listofshocks)
         innovations(1,:) = innovations(1,:) + transpose(C(:,j));
     end
     if options_.linear
-        ysim__1 = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic);
+        [ysim__1, errorflag] = simul_backward_linear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic);
     else
-        ysim__1 = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic);
+        [ysim__1, errorflag] = simul_backward_nonlinear_model_(initialcondition, periods, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic);
+    end
+    if errorflag
+        warning('Simulation failed. Cannot compute IRF for %s.', listofshocks{i})
+        continue
     end
     % Transform the endogenous variables
     if notransform
@@ -230,4 +238,4 @@ end
 if nargout>1
     baseline = dseries(transpose(endo_simul__0), initialcondition.init, endonames(1:M_.orig_endo_nbr), DynareModel.endo_names_tex(1:M_.orig_endo_nbr));
     baseline = merge(baseline, innovationbaseline);
-end
\ No newline at end of file
+end
diff --git a/matlab/backward/simul_backward_linear_model.m b/matlab/backward/simul_backward_linear_model.m
index ea7dc2758a8d3a89fabee5ffd537a80d8aea36a2..2c9252513ddc8da9cf72eed5b34cbc6ec7c8962c 100644
--- a/matlab/backward/simul_backward_linear_model.m
+++ b/matlab/backward/simul_backward_linear_model.m
@@ -1,4 +1,4 @@
-function simulations = simul_backward_linear_model(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations)
+function [simulations, errorflag] = simul_backward_linear_model(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations)
 
 % Simulates a stochastic linear backward looking model.
 %
@@ -12,6 +12,7 @@ function simulations = simul_backward_linear_model(initialconditions, samplesize
 %
 % OUTPUTS
 % - DynareOutput        [struct]      Dynare's oo_ global structure.
+% - errorflag           [logical]     scalar, equal to false iff the simulation did not fail.
 %
 % REMARKS
 % [1] The innovations used for the simulation are saved in DynareOutput.exo_simul, and the resulting paths for the endogenous
@@ -21,7 +22,7 @@ function simulations = simul_backward_linear_model(initialconditions, samplesize
 % [3] If the first input argument is empty, the endogenous variables are initialized with 0, or if available with the informations
 %     provided thrtough the histval block.
 
-% Copyright © 2012-2022 Dynare Team
+% Copyright © 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -53,6 +54,6 @@ end
 [initialconditions, samplesize, innovations, DynareOptions, DynareModel, DynareOutput, endonames, exonames, nx, ny1, iy1, jdx, model_dynamic] = ...
     simul_backward_model_init(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations);
 
-[ysim, xsim] = simul_backward_linear_model_(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic);
+[ysim, xsim, errorflag] = simul_backward_linear_model_(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic);
 
-simulations = [dseries(ysim', initialconditions.init, endonames(1:DynareModel.orig_endo_nbr)), dseries(xsim, initialconditions.init, exonames)];
\ No newline at end of file
+simulations = [dseries(ysim', initialconditions.init, endonames(1:DynareModel.orig_endo_nbr)), dseries(xsim, initialconditions.init, exonames)];
diff --git a/matlab/backward/simul_backward_linear_model_.m b/matlab/backward/simul_backward_linear_model_.m
index 12f51254ab9b03a590627ae36a0543ace673e284..3f71ac1eb92ec293f2d23eada3ac92f903d467c6 100644
--- a/matlab/backward/simul_backward_linear_model_.m
+++ b/matlab/backward/simul_backward_linear_model_.m
@@ -1,4 +1,4 @@
-function [ysim, xsim] = simul_backward_linear_model_(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic)
+function [ysim, xsim, errorflag] = simul_backward_linear_model_(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations, nx, ny1, iy1, jdx, model_dynamic)
 
 % Simulates a stochastic linear backward looking model.
 %
@@ -12,6 +12,7 @@ function [ysim, xsim] = simul_backward_linear_model_(initialconditions, samplesi
 %
 % OUTPUTS
 % - DynareOutput        [struct]      Dynare's oo_ global structure.
+% - errorflag           [logical]     scalar, equal to false iff the simulation did not fail.
 %
 % REMARKS
 % [1] The innovations used for the simulation are saved in DynareOutput.exo_simul, and the resulting paths for the endogenous
@@ -21,7 +22,7 @@ function [ysim, xsim] = simul_backward_linear_model_(initialconditions, samplesi
 % [3] If the first input argument is empty, the endogenous variables are initialized with 0, or if available with the informations
 %     provided thrtough the histval block.
 
-% Copyright © 2017-2022 Dynare Team
+% Copyright © 2017-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -38,6 +39,8 @@ function [ysim, xsim] = simul_backward_linear_model_(initialconditions, samplesi
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
+errorflag = false;
+
 if ~isempty(innovations)
     DynareOutput.exo_simul(initialconditions.nobs+(1:samplesize),:) = innovations;
 end
@@ -48,7 +51,15 @@ end
                              DynareModel.params, ...
                              DynareOutput.steady_state, DynareModel.orig_maximum_lag+1);
 
-A0inv = inv(jacob(:,jdx));
+try
+    A0inv = inv(jacob(:,jdx));
+catch
+    errorflag = true;
+    ysim = [];
+    xsim = [];
+    return
+end
+
 A1 = jacob(:,nonzeros(DynareModel.lead_lag_incidence(1,:)));
 B = jacob(:,end-nx+1:end);
 
@@ -58,4 +69,4 @@ for it = initialconditions.nobs+(1:samplesize)
 end
 
 ysim = DynareOutput.endo_simul(1:DynareModel.orig_endo_nbr,:);
-xsim = DynareOutput.exo_simul;
\ No newline at end of file
+xsim = DynareOutput.exo_simul;
diff --git a/matlab/backward/simul_backward_model.m b/matlab/backward/simul_backward_model.m
index eb4f13088842dc7baca20d4439a4a54359e77ea2..0cf99c18e3082f108929a3df91a6c772c6584f3f 100644
--- a/matlab/backward/simul_backward_model.m
+++ b/matlab/backward/simul_backward_model.m
@@ -1,4 +1,4 @@
-function simulation = simul_backward_model(initialconditions, samplesize, innovations)
+function [simulation, errorflag] = simul_backward_model(initialconditions, samplesize, innovations)
 % function simulation = simul_backward_model(initialconditions, samplesize, innovations)
 % Simulates a stochastic backward looking model (with arbitrary precision).
 %
@@ -9,6 +9,7 @@ function simulation = simul_backward_model(initialconditions, samplesize, innova
 %
 % OUTPUTS
 % - simulation          [dseries]     Simulated endogenous and exogenous variables.
+% - errorflag           [logical]     scalar, equal to false iff the simulation did not fail.
 %
 % REMARKS
 % [1] The innovations used for the simulation are saved in DynareOutput.exo_simul, and the resulting paths for the endogenous
@@ -18,7 +19,7 @@ function simulation = simul_backward_model(initialconditions, samplesize, innova
 % [3] If the first input argument is empty, the endogenous variables are initialized with 0, or if available with the informations
 %     provided through the histval block.
 
-% Copyright © 2012-2022 Dynare Team
+% Copyright © 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -93,7 +94,7 @@ else
 end
 
 if options_.linear
-    simulation = simul_backward_linear_model(initialconditions, samplesize, options_, M_, oo_, Innovations);
+    [simulation, errorflag] = simul_backward_linear_model(initialconditions, samplesize, options_, M_, oo_, Innovations);
 else
-    simulation = simul_backward_nonlinear_model(initialconditions, samplesize, options_, M_, oo_, Innovations);
+    [simulation, errorflag] = simul_backward_nonlinear_model(initialconditions, samplesize, options_, M_, oo_, Innovations);
 end
diff --git a/matlab/backward/simul_backward_nonlinear_model.m b/matlab/backward/simul_backward_nonlinear_model.m
index d4df3c39dfcfca9c7b3f4c3f9a2e2f2840b85409..d1512c7b7e53140eeab447c8996d1f38f496e9fa 100644
--- a/matlab/backward/simul_backward_nonlinear_model.m
+++ b/matlab/backward/simul_backward_nonlinear_model.m
@@ -1,4 +1,4 @@
-function simulations = simul_backward_nonlinear_model(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations)
+function [simulations, errorflag] = simul_backward_nonlinear_model(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations)
 % function simulations = simul_backward_nonlinear_model(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations)
 % Simulates a stochastic non linear backward looking model with arbitrary precision (a deterministic solver is used).
 %
@@ -12,6 +12,7 @@ function simulations = simul_backward_nonlinear_model(initialconditions, samples
 %
 % OUTPUTS
 % - simulation          [dseries]     Simulated endogenous and exogenous variables.
+% - errorflag           [logical]     scalar, equal to false iff the simulation did not fail.
 %
 % REMARKS
 % [1] The innovations used for the simulation are saved in DynareOutput.exo_simul, and the resulting paths for the endogenous
@@ -21,7 +22,7 @@ function simulations = simul_backward_nonlinear_model(initialconditions, samples
 % [3] If the first input argument is empty, the endogenous variables are initialized with 0, or if available with the informations
 %     provided thrtough the histval block.
 
-% Copyright (©) 2012-2022 Dynare Team
+% Copyright (©) 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -53,6 +54,6 @@ end
 [initialconditions, samplesize, innovations, DynareOptions, DynareModel, DynareOutput, endonames, exonames, ~, ~, iy1, ~, model_dynamic] = ...
     simul_backward_model_init(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations);
 
-[ysim, xsim] = simul_backward_nonlinear_model_(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic);
+[ysim, xsim, errorflag] = simul_backward_nonlinear_model_(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic);
 
-simulations = [dseries(ysim', initialconditions.init, endonames(1:DynareModel.orig_endo_nbr)), dseries(xsim, initialconditions.init, exonames)];
\ No newline at end of file
+simulations = [dseries(ysim', initialconditions.init, endonames(1:DynareModel.orig_endo_nbr)), dseries(xsim, initialconditions.init, exonames)];
diff --git a/matlab/backward/simul_backward_nonlinear_model_.m b/matlab/backward/simul_backward_nonlinear_model_.m
index 024b8f1730ec8885f33c39bec30d0a4c6bce3936..146c878fce4f265d92438b721b1da57c6ccf83a6 100644
--- a/matlab/backward/simul_backward_nonlinear_model_.m
+++ b/matlab/backward/simul_backward_nonlinear_model_.m
@@ -1,4 +1,4 @@
-function [ysim, xsim] = simul_backward_nonlinear_model_(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic)
+function [ysim, xsim, errorflag] = simul_backward_nonlinear_model_(initialconditions, samplesize, DynareOptions, DynareModel, DynareOutput, innovations, iy1, model_dynamic)
 
 % Simulates a stochastic non linear backward looking model with arbitrary precision (a deterministic solver is used).
 %
@@ -12,6 +12,7 @@ function [ysim, xsim] = simul_backward_nonlinear_model_(initialconditions, sampl
 %
 % OUTPUTS
 % - DynareOutput        [struct]      Dynare's oo_ global structure.
+% - errorflag           [logical]     scalar, equal to false iff the simulation did not fail.
 %
 % REMARKS
 % [1] The innovations used for the simulation are saved in DynareOutput.exo_simul, and the resulting paths for the endogenous
@@ -21,7 +22,7 @@ function [ysim, xsim] = simul_backward_nonlinear_model_(initialconditions, sampl
 % [3] If the first input argument is empty, the endogenous variables are initialized with 0, or if available with the informations
 %     provided thrtough the histval block.
 
-% Copyright © 2017-2022 Dynare Team
+% Copyright © 2017-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -100,6 +101,7 @@ for it = initialconditions.nobs+(1:samplesize)
             end
         end
     catch Error
+        errorflag = true;
         DynareOutput.endo_simul = DynareOutput.endo_simul(:, 1:it-1);
         dprintf('Newton failed on iteration i = %s.', num2str(it-initialconditions.nobs));
         ytm = DynareOutput.endo_simul(:,end);
diff --git a/tests/ecb/backward-models/solow_forecast.mod b/tests/ecb/backward-models/solow_forecast.mod
index 6832141423c0d9839e6a4588c630f00358a13e08..78d766ec9d03446c7e5ad20f7ae906a0ab3c026f 100644
--- a/tests/ecb/backward-models/solow_forecast.mod
+++ b/tests/ecb/backward-models/solow_forecast.mod
@@ -51,7 +51,7 @@ histval;
     PhysicalCapitalStock(0) = 1;
 end;
 
-simulation__ = simul_backward_nonlinear_model([], 10, options_, M_, oo_);
+[simulation__, errorflag] = simul_backward_nonlinear_model([], 10, options_, M_, oo_);
 
 initialcondition = dseries(simulation__.data(10,1:M_.endo_nbr), 2017Q1, M_.endo_names);
 
@@ -73,11 +73,11 @@ listofvariables = {'Efficiency', 'Output'};
 forecasts = backward_model_forecast(initialcondition, listofvariables, 16, true);
 
 /* REMARKS
-** 
+**
 ** - The third argument is the number of periods for the forecast. Default is 8.
 ** - The last argument is for computing (or not) the confidence bands. If false (default), only point forecast is produced.
 ** - forecasts is a structure, each field is a dseries object for the point forecast (ie forecasts without innovations in the future),
-** the mean forecast, the median forecast, the standard deviation of the predictive distribution and the lower/upper bounds of the 
+** the mean forecast, the median forecast, the standard deviation of the predictive distribution and the lower/upper bounds of the
 ** interval containing 95% of the predictive distribution.
 */
 
@@ -89,7 +89,7 @@ plot(forecasts.ub.Output,'-r')
 hold off
 
 /* REMARKS
-** 
-** In this model there is no steady state (only a stable balanced growth paths), this explains the 
+**
+** In this model there is no steady state (only a stable balanced growth paths), this explains the
 ** shape of the forecast for Output.
-*/
\ No newline at end of file
+*/