diff --git a/license.txt b/license.txt
index d5488b35a..b00407108 100644
--- a/license.txt
+++ b/license.txt
@@ -119,7 +119,7 @@ License: GPL-3+
 
 Files: matlab/trust_region.m
 Copyright: 2008-2012 VZLU Prague, a.s.
-           2014-2019 Dynare Team
+           2014-2021 Dynare Team
 License: GPL-3+
 
 Files: matlab/one_sided_hp_filter.m
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 85d15d507..794554f3b 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -232,9 +232,10 @@ elseif options.solve_algo==1
                             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, ...
-                             maxit, options.debug, arguments{:});
+    [x, errorflag] = 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{:});
 elseif ismember(options.solve_algo, [2, 12, 4])
     if ismember(options.solve_algo, [2, 12])
         solver = @solve1;
@@ -312,8 +313,9 @@ elseif ismember(options.solve_algo, [2, 12, 4])
         end
         [x, errorflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
                                 options.gstep, ...
-                                tolf, options.solve_tolx, ...
-                                maxit, options.debug, arguments{:});
+                                tolf, options.solve_tolx, maxit, ...
+                                options.trust_region_initial_step_bound_factor, ...
+                                options.debug, arguments{:});
         fre = true;
         if errorflag
             return
@@ -323,8 +325,9 @@ elseif ismember(options.solve_algo, [2, 12, 4])
     if max(abs(fvec))>tolf
         disp('Call solver on the full nonlinear problem.')
         [x, errorflag] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
-                                options.gstep, tolf, options.solve_tolx, ...
-                                maxit, options.debug, arguments{:});
+                                options.gstep, tolf, options.solve_tolx, maxit, ...
+                                options.trust_region_initial_step_bound_factor, ...
+                                options.debug, arguments{:});
     end
 elseif options.solve_algo==3
     if jacobian_flag
diff --git a/matlab/perfect-foresight-models/sim1_purely_forward.m b/matlab/perfect-foresight-models/sim1_purely_forward.m
index 130989fe1..ff893bd4e 100644
--- a/matlab/perfect-foresight-models/sim1_purely_forward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_forward.m
@@ -32,9 +32,10 @@ info.status = 1;
 for it = options.periods:-1:1
     yf = endogenousvariables(:,it+1); % Values at next period, also used as guess value for current period
     yf1 = yf(iyf);
+    % TODO Call dynare_solve instead.
     [tmp, check] = solve1(dynamicmodel, [yf; yf1], 1:M.endo_nbr, 1:M.endo_nbr, ...
                           1, options.gstep, options.dynatol.f, ...
-                          options.dynatol.x, options.simul.maxit, ...
+                          options.dynatol.x, options.simul.maxit, [], ...
                           options.debug, exogenousvariables, M.params, steadystate, ...
                           it+M.maximum_lag);
     if check
diff --git a/matlab/solve1.m b/matlab/solve1.m
index 73dca70b2..8a1421dd6 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -1,4 +1,4 @@
-function [x,check] = solve1(func,x,j1,j2,jacobian_flag,gstep,tolf,tolx,maxit,debug,varargin)
+function [x,check] = solve1(func,x,j1,j2,jacobian_flag,gstep,tolf,tolx,maxit,fake,debug,varargin)
 % Solves systems of non linear equations of several variables
 %
 % INPUTS
@@ -13,6 +13,7 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,gstep,tolf,tolx,maxit,deb
 %    tolf             tolerance for residuals
 %    tolx             tolerance for solution variation
 %    maxit            maximum number of iterations
+%    fake             unused argument (compatibity with trust_region).
 %    debug            debug flag
 %    varargin:        list of extra arguments to the function
 %
@@ -23,7 +24,7 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,gstep,tolf,tolx,maxit,deb
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2001-2017 Dynare Team
+% Copyright © 2001-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
diff --git a/matlab/trust_region.m b/matlab/trust_region.m
index d5446803c..ab5d6a908 100644
--- a/matlab/trust_region.m
+++ b/matlab/trust_region.m
@@ -1,4 +1,4 @@
-function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,debug,varargin)
+function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,factor,debug,varargin)
 % Solves systems of non linear equations of several variables, using a
 % trust-region method.
 %
@@ -14,6 +14,7 @@ function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tol
 %    tolf             tolerance for residuals
 %    tolx             tolerance for solution variation
 %    maxiter          maximum number of iterations
+%    factor           real scalar, determines the initial step bound
 %    debug            debug flag
 %    varargin:        list of arguments following bad_cond_flag
 %
@@ -101,7 +102,7 @@ while (niter < maxiter && ~info)
     if (niter == 1)
         xn = norm (dg .* x(j2));
         % FIXME: something better?
-        delta = max (xn, 1);
+        delta = max (xn, 1)*factor;
     end
 
     % Get trust-region model (dogleg) minimizer.