diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index aa9d5fd3aab87989bba1512eb6652d198a648913..1beb8f6311c6600d802a3298c91190d8b1587d89 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -2954,20 +2954,17 @@ Finding the steady state with Dynare nonlinear solver
``6``
Newton algorithm with a sparse LU solver at each
- iteration (requires ``bytecode`` and/or ``block``
- option, see :ref:`model-decl`).
+ iteration.
``7``
Newton algorithm with a Generalized Minimal Residual
- (GMRES) solver at each iteration (requires ``bytecode``
- and/or ``block`` option, see :ref:`model-decl`).
+ (GMRES) solver at each iteration.
``8``
Newton algorithm with a Stabilized Bi-Conjugate
- Gradient (BiCGStab) solver at each iteration (requires
- bytecode and/or block option, see :ref:`model-decl`).
+ Gradient (BiCGStab) solver at each iteration.
``9``
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 292ef918c5294f20a65d9ef49eb65750b12ab043..6e78515fa22fb850e4190b9ac1dcea9ebda0bf46 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,exo_ss,M_,options_,ste
% SPECIAL REQUIREMENTS
% none
-% Copyright © 2001-2023 Dynare Team
+% Copyright © 2001-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -43,11 +43,6 @@ if options_.solve_algo < 0 || options_.solve_algo > 14
error('STEADY: solve_algo must be between 0 and 14')
end
-if ~options_.bytecode && ~options_.block && options_.solve_algo > 4 && ...
- options_.solve_algo < 9
- error('STEADY: you can''t use solve_algo = {5,6,7,8} without block nor bytecode options_')
-end
-
if ~options_.bytecode && options_.block && options_.solve_algo == 5
error('STEADY: you can''t use solve_algo = 5 without bytecode option')
end
diff --git a/matlab/optimization/dynare_solve.m b/matlab/optimization/dynare_solve.m
index 1278055b215e3029a77853b6542f776dc3047538..3779aae7c091a15a2c0e63895bc4e89d8b377700 100644
--- a/matlab/optimization/dynare_solve.m
+++ b/matlab/optimization/dynare_solve.m
@@ -293,6 +293,9 @@ elseif options_.solve_algo==3
errorflag = true;
end
[fvec, fjac] = feval(f, x, varargin{:});
+elseif ismember(options_.solve_algo, [6 7 8])
+ [x, errorflag, errorcode] = newton_solve(f, x, jacobian_flag, options_.gstep, tolf, tolx, maxit, options_.solve_algo, varargin{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
elseif options_.solve_algo==10
% LMMCP
olmmcp = options_.lmmcp;
@@ -337,5 +340,5 @@ elseif ismember(options_.solve_algo, [13, 14])
options_.debug, varargin{:});
[fvec, fjac] = feval(f, x, varargin{:});
else
- error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,7,8,9,10,11,12,13,14]')
+ error('DYNARE_SOLVE: option solve_algo must be one of [0:4,6:14]')
end
diff --git a/matlab/optimization/newton_solve.m b/matlab/optimization/newton_solve.m
new file mode 100644
index 0000000000000000000000000000000000000000..de1d2e48286ebd95973d632d1d7d3873e7f9a3bb
--- /dev/null
+++ b/matlab/optimization/newton_solve.m
@@ -0,0 +1,131 @@
+function [x, errorflag, errorcode] = newton_solve(func, x, jacobian_flag, gstep, tolf, tolx, maxit, solve_algo, varargin)
+
+% Solves systems of non linear equations of several variables using a Newton solver, with three
+% variants for the inner linear solver:
+% solve_algo = 6: use a sparse LU
+% solve_algo = 7: use GMRES
+% solve_algo = 8: use BiCGStab
+%
+% INPUTS
+% func: name of the function to be solved
+% x: guess values
+% jacobian_flag=true: jacobian given by the 'func' function
+% jacobian_flag=false: jacobian obtained numerically
+% gstep increment multiplier in numerical derivative
+% computation
+% tolf tolerance for residuals
+% tolx tolerance for solution variation
+% maxit maximum number of iterations
+% solve_algo from options_
+% varargin: list of extra arguments to the function
+%
+% OUTPUTS
+% x: results
+% errorflag=true: the model can not be solved
+
+% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
+
+nn = length(x);
+
+errorflag = false;
+
+% Declaration needed for sharing with nested function
+fvec = NaN(nn);
+fjac = NaN(nn, nn);
+
+function update_fvec_fjac
+ if jacobian_flag
+ [fvec, fjac] = feval(func, x, varargin{:});
+ else
+ fvec = feval(func, x, varargin{:});
+ dh = max(abs(x),gstep(1)*ones(nn,1))*eps^(1/3);
+ for j = 1:nn
+ xdh = x;
+ xdh(j) = xdh(j)+dh(j);
+ t = feval(func, xdh, varargin{:});
+ fjac(:, j) = (t - fvec) ./ dh(j);
+ end
+ end
+end
+
+update_fvec_fjac;
+
+idInf = isinf(fvec);
+idNan = isnan(fvec);
+idCpx = ~isreal(fvec);
+
+if any(idInf)
+ disp('newton_solve: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number:')
+ disp(idInf')
+ errorcode = 0;
+ errorflag = true;
+ return
+end
+
+if any(idNan)
+ disp('newton_solve: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a nan:')
+ disp(idNan')
+ errorcode = 0;
+ errorflag = true;
+ return
+end
+
+if any(idNan)
+ disp('newton_solve: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a complex number:')
+ disp(idCpx')
+ errorcode = 0;
+ errorflag = true;
+ return
+end
+
+if max(abs(fvec)) < tolf
+ % Initial guess is a solution
+ errorcode = -1;
+ return
+end
+
+for it = 1:maxit
+ if solve_algo == 6
+ p = fjac\fvec;
+ else
+ ilu_setup.type = 'ilutp';
+ ilu_setup.droptol = 1e-10;
+ [L1, U1] = ilu(fjac, ilu_setup);
+ if solve_algo == 7
+ p = gmres(fjac, fvec, [], [], [], L1, U1);
+ elseif solve_algo == 8
+ p = bicgstab(fjac, fvec, [], [], L1, U1);
+ else
+ error('newton_solve: invalid value for solve_algo')
+ end
+ end
+
+ x = x - p;
+
+ update_fvec_fjac;
+
+ if max(abs(fvec)) < tolf
+ errorcode = 1;
+ return
+ end
+end
+
+errorflag = true;
+errorcode = 2;
+
+end
diff --git a/tests/run_block_bytecode_tests.m b/tests/run_block_bytecode_tests.m
index e9aac919470257f38ae3f8b3fbad02e8272cc463..0c95b794ab539939b7c58eefbb4d0d51b730c674 100644
--- a/tests/run_block_bytecode_tests.m
+++ b/tests/run_block_bytecode_tests.m
@@ -41,7 +41,7 @@ for blockFlag = 0:1
default_solve_algo = 2;
default_stack_solve_algo = 0;
if ~blockFlag && storageFlag ~= 2
- solve_algos = [1:4 9];
+ solve_algos = [1:4 6:9];
stack_solve_algos = [0:3 6];
elseif blockFlag && storageFlag ~= 2
solve_algos = [1:4 6:9];