diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index b27e3d91eb7109dd99c97e496af5931e38860858..57b8288a273bbfea8e0a2b81a9a96134a77d89ed 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -15,7 +15,7 @@ function [x, errorflag, fvec, fjac] = dynare_solve(f, x, options, varargin)
 % - 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).
 
-% Copyright © 2001-2020 Dynare Team
+% Copyright © 2001-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -225,8 +225,8 @@ if options.solve_algo == 0
     end
 elseif options.solve_algo==1
     [x, errorflag] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, ...
-                       tolf, tolx, ...
-                       maxit, options.debug, arguments{:});
+                            tolf, tolx, ...
+                            maxit, options.debug, arguments{:});
 elseif options.solve_algo==9
     [x, errorflag] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, ...
                              tolf, tolx, ...
diff --git a/matlab/perfect-foresight-models/sim1_purely_backward.m b/matlab/perfect-foresight-models/sim1_purely_backward.m
index 689d4f310d8e490eaca3bf48d5e195a92c1cd9db..daee720cea7e62dafb9d8b869ca7b265218f8c72 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 (C) 2012-2017 Dynare Team
+% Copyright © 2012-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -29,24 +29,31 @@ else
     iyb = [];
 end
 
-
 if ny0 ~= M.endo_nbr
     error('All endogenous variables must appear at the current period!')
 end
 
+if ismember(options.solve_algo, [12,14]) && ~M.possible_to_use_solve_algo_12_14
+    error(M.message_solve_algo_12_14)
+end
+
 dynamicmodel = str2func([M.fname,'.dynamic']);
+dynamicmodel_s = str2func('dynamic_backward_model_for_simulation');
 
-info.status = 1;
+info.status = true;
 
 for it = M.maximum_lag + (1:options.periods)
-    yb = endogenousvariables(:,it-1);        % Values at previous period, also used as guess value for current period
-    yb1 = yb(iyb);
-    [tmp, check] = solve1(dynamicmodel, [yb1; yb], 1:M.endo_nbr, nyb+1:nyb+M.endo_nbr, ...
-                          1, options.gstep, options.dynatol.f, options.dynatol.x, ...
-                          options.simul.maxit, options.debug, exogenousvariables, ...
-                          M.params, steadystate, it);
+    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);
+    else
+        [tmp, check] = dynare_solve(dynamicmodel_s, y, options, ...
+                                    dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
+    end
     if check
-        info.status = 0;
+        info.status = false;
     end
-    endogenousvariables(:,it) = tmp(nyb+1:nyb+M.endo_nbr);
+    endogenousvariables(:,it) = tmp;
 end
\ No newline at end of file
diff --git a/tests/simul/Solow_no_varexo.mod b/tests/simul/Solow_no_varexo.mod
index 542fd32584418e166e2abe9d8aacc3b988c3365d..f790061d0ec6b22d5018d7a77cea05231f36f096 100644
--- a/tests/simul/Solow_no_varexo.mod
+++ b/tests/simul/Solow_no_varexo.mod
@@ -61,5 +61,5 @@ end;
 resid(1);
 steady;
 
-simul(periods=200);
+simul(periods=200, solve_algo=4);
 rplot k;