diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 929079add12b208366a19ca4451c5d887a9cb333..2f4c7bb18c00baeff986052cbc6ac585b38adeca 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -328,6 +328,17 @@ options_.simul.endval_steady = false;
 options_.simul.first_simulation_period = dates();
 options_.simul.last_simulation_period = dates();
 
+% For stack_solve_algo={2,3}
+options_.simul.preconditioner = 'iterstack';
+options_.simul.itertol = 1e-12;
+options_.simul.itermaxit = 500;
+options_.simul.gmres_restart = 100;
+options_.simul.iterstack.maxlu = 20000; % Similar to TROLL’s ITERSTACK
+options_.simul.iterstack.relu = 0.5; % Similar to TROLL’s ITERSTACK
+options_.simul.ilu.type = 'ilutp';
+options_.simul.ilu.droptol = 1e-12;
+options_.simul.ilu.udiag = true;
+
 options_.simul.homotopy_max_completion_share = 1;
 options_.simul.homotopy_min_step_size = 1e-3;
 options_.simul.homotopy_step_size_increase_success_count = 3;
diff --git a/matlab/optimization/solve_one_boundary.m b/matlab/optimization/solve_one_boundary.m
index e9177db7d02acb3d508cc5eb9a7ef8d566b17c3d..3c4ca9a6cc98cc1026514fa4142df25552c13a70 100644
--- a/matlab/optimization/solve_one_boundary.m
+++ b/matlab/optimization/solve_one_boundary.m
@@ -40,7 +40,7 @@ function [y, T, success, max_res, iter] = solve_one_boundary(fh, y, x, params, s
 % ALGORITHM
 %   Newton with LU or GMRES or BicGstab for dynamic block
 
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -59,11 +59,8 @@ function [y, T, success, max_res, iter] = solve_one_boundary(fh, y, x, params, s
 
 Blck_size=size(y_index_eq,2);
 correcting_factor=0.01;
-ilu_setup.type='crout';
-ilu_setup.droptol=1e-10;
 max_resa=1e100;
 lambda = 1; % Length of Newton step
-reduced = 0;
 if is_forward
     incr = 1;
     start = y_kmin+1;
@@ -141,7 +138,6 @@ for it_=start:incr:finish
                         end
                     elseif lambda>1e-8
                         lambda=lambda/2;
-                        reduced = 1;
                         if verbose
                             disp(['reducing the path length: lambda=' num2str(lambda,'%f')])
                         end
@@ -192,7 +188,7 @@ for it_=start:incr:finish
                                     M_.block_structure.block(Block_Num).g1_sparse_rowval, ...
                                     M_.block_structure.block(Block_Num).g1_sparse_colval, ...
                                     M_.block_structure.block(Block_Num).g1_sparse_colptr, T(:, it_));
-            elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0 || stack_solve_algo==6)) || (~is_dynamic && options_.solve_algo==6)
+            elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0 || stack_solve_algo==6 || (ismember(stack_solve_algo, [2, 3]) && strcmp(options_.simul.preconditioner, 'iterstack') && size(g1, 1) < options_.simul.iterstack.maxlu))) || (~is_dynamic && options_.solve_algo==6) % Fallback to LU if block too small for iterstack
                 if verbose && ~is_dynamic
                     disp('steady: Sparse LU ')
                 end
@@ -203,64 +199,37 @@ for it_=start:incr:finish
                 else
                     y(y_index_eq) = ya;
                 end
-            elseif (stack_solve_algo==2 && is_dynamic) || (options_.solve_algo==7 && ~is_dynamic)
-                flag1=1;
+            elseif (is_dynamic && ismember(stack_solve_algo, [2 3])) || (~is_dynamic && ismember(options_.solve_algo, [7 8]))
                 if verbose && ~is_dynamic
-                    disp('steady: GMRES ')
-                end
-                while flag1>0
-                    [L1, U1]=ilu(g1,ilu_setup);
-                    [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
-                    if  flag1>0 || reduced
-                        if verbose
-                            if flag1==1
-                                disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')])
-                            elseif(flag1==2)
-                                disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')])
-                            elseif(flag1==3)
-                                disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')])
-                            end
-                        end
-                        ilu_setup.droptol = ilu_setup.droptol/10;
-                        reduced = 0;
+                    if options_.solve_algo == 7
+                        disp('steady: GMRES ')
                     else
-                        ya = ya + lambda*dx;
-                        if is_dynamic
-                            y(y_index_eq, it_) = ya;
-                        else
-                            y(y_index_eq) = ya';
-                        end
+                        disp('steady: BiCGStab')
                     end
                 end
-            elseif (stack_solve_algo==3 && is_dynamic) || (options_.solve_algo==8 && ~is_dynamic)
-                flag1=1;
-                if verbose && ~is_dynamic
-                    disp('steady: BiCGStab')
+                if strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0
+                    [L, U, P, Q] = lu(g1);
+                    zc = L\(P*r);
+                    zb = U\zc;
+                elseif strcmp(options_.simul.preconditioner, 'iterstack')
+                    [L, U, P, Q] = iterstack_preconditioner(g1, options_);
+                elseif strcmp(options_.simul.preconditioner, 'ilu')
+                    [L, U, P] = ilu(g1, options_.simul.ilu);
+                    Q = speye(size(g1));
                 end
-                while flag1>0
-                    [L1, U1]=ilu(g1,ilu_setup);
-                    [dx,flag1] = bicgstab(g1,-r,1e-6,Blck_size,L1,U1);
-                    if flag1>0 || reduced
-                        if verbose
-                            if(flag1==1)
-                                disp(['Error in simul: No convergence inside BiCGStab after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')])
-                            elseif(flag1==2)
-                                disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')])
-                            elseif(flag1==3)
-                                disp(['Error in simul: BiCGStab stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')])
-                            end
-                        end
-                        ilu_setup.droptol = ilu_setup.droptol/10;
-                        reduced = 0;
+                if ~(strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0)
+                    if (stack_solve_algo==2 && is_dynamic) || (options_.solve_algo==7 && ~is_dynamic)
+                        zb = gmres(P*g1a*Q, P*r, options_.simul.gmres_restart, options_.simul.itertol, options_.simul.itermaxit, L, U);
                     else
-                        ya = ya + lambda*dx;
-                        if is_dynamic
-                            y(y_index_eq, it_) = ya;
-                        else
-                            y(y_index_eq) = ya';
-                        end
+                        zb = bicgstab(P*g1a*Q, P*r, options_.simul.itertol, options_.simul.itermaxit, L, U);
                     end
                 end
+                ya = ya - lambda*Q*zb;
+                if is_dynamic
+                    y(y_index_eq, it_) = ya;
+                else
+                    y(y_index_eq) = ya';
+                end
             else
                 if is_dynamic
                     error(['options_.stack_solve_algo = ' num2str(stack_solve_algo) ' not implemented'])
diff --git a/matlab/perfect-foresight-models/iterstack_preconditioner.m b/matlab/perfect-foresight-models/iterstack_preconditioner.m
new file mode 100644
index 0000000000000000000000000000000000000000..ac3205c70906b806c63c260fb5385963a1ad0cf0
--- /dev/null
+++ b/matlab/perfect-foresight-models/iterstack_preconditioner.m
@@ -0,0 +1,41 @@
+function [L, U, P, Q] = iterstack_preconditioner(A, options_)
+% Copyright © 2025 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/>.
+
+if size(A, 1) < options_.simul.iterstack.maxlu
+    error('iterstack preconditionner: too small problem') % TODO: better recovery
+end
+ny = size(A, 1) / options_.periods; % TODO: should ideally be passed as argument
+if ny > options_.simul.iterstack.maxlu
+    error('iterstack preconditionner: too large model; either increase maxlu or switch to another algorithm') % TODO: improve message
+end
+lusize = floor(options_.simul.iterstack.maxlu / ny) * ny;
+luidx = floor((floor(size(A, 1) / lusize)-1) * options_.simul.iterstack.relu) * lusize + (1:lusize);
+[L1,U1,P1,Q1] = lu(A(luidx,luidx));
+eyek = speye(floor(size(A, 1) / lusize));
+L = kron(eyek, L1);
+U = kron(eyek, U1);
+P = kron(eyek, P1);
+Q = kron(eyek, Q1);
+r = rem(size(A, 1), lusize);
+if r > 0 % Compute additional smaller LU for remainder, if any
+    [L2,U2,P2,Q2] = lu(A(end-r+1:end,end-r+1:end));
+    L = blkdiag(L, L2);
+    U = blkdiag(U, U2);
+    P = blkdiag(P, P2);
+    Q = blkdiag(Q, Q2);
+end
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 974cda2d5ccdabf650ba0a9384e0b187b6383db6..42be5ec776bfb96e30b48867f9e698190e7063d0 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -148,7 +148,7 @@ while ~(converged || iter > options_.simul.maxit)
         if options_.simul.robust_lin_solve
             dy = -lin_solve_robust(A, res, verbose, options_);
         else
-            dy = -lin_solve(A, res, verbose, options_);
+            dy = -lin_solve(A, res, verbose, options_, iter == 1);
         end
         if any(isnan(dy)) || any(isinf(dy))
             if verbose
@@ -217,7 +217,10 @@ if verbose
     skipline();
 end
 
-function x = lin_solve(A, b, verbose, options_)
+function x = lin_solve(A, b, verbose, options_, first_iter)
+
+persistent L U P Q
+
 if norm(b) < sqrt(eps) % then x = 0 is a solution
     x = 0;
     return
@@ -226,16 +229,26 @@ end
 if options_.stack_solve_algo == 0
     x = A\b;
 else
-    ilu_setup.type = 'ilutp';
-    ilu_setup.droptol = 1e-10;
-    [L1, U1] = ilu(A, ilu_setup);
-    if options_.stack_solve_algo == 2
-        x = gmres(A, b, [], [], [], L1, U1);
-    elseif options_.stack_solve_algo == 3
-        x = bicgstab(A, b, [], [], L1, U1);
-    else
-        error('sim1: invalid value for options_.stack_solve_algo')
+    if strcmp(options_.simul.preconditioner, 'umfiter') && first_iter
+        [L, U, P, Q] = lu(A);
+        z = L\(P*b);
+        y = U\z;
+    elseif strcmp(options_.simul.preconditioner, 'iterstack')
+        [L, U, P, Q] = iterstack_preconditioner(A, options_);
+    elseif strcmp(options_.simul.preconditioner, 'ilu')
+        [L, U, P] = ilu(A, options_.simul.ilu);
+        Q = speye(size(A));
+    end
+    if ~(strcmp(options_.simul.preconditioner, 'umfiter') && first_iter)
+        if options_.stack_solve_algo == 2
+            y = gmres(P*A*Q, P*b, options_.simul.gmres_restart, options_.simul.itertol, options_.simul.itermaxit, L, U);
+        elseif options_.stack_solve_algo == 3
+            y = bicgstab(P*A*Q, P*b, options_.simul.itertol, options_.simul.itermaxit, L, U);
+        else
+            error('sim1: invalid value for options_.stack_solve_algo')
+        end
     end
+    x = Q*y;
 end
 x(~isfinite(x)) = 0;
 relres = norm(b - A*x) / norm(b);
diff --git a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
index 14db6dee693ee24108b535248d5c5b5147395361..16840c89f426b893d1640c8b65e5c7e85f87785f 100644
--- a/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
+++ b/matlab/perfect-foresight-models/solve_block_decomposed_problem.m
@@ -40,9 +40,9 @@ switch options_.stack_solve_algo
     case {1,6}
         mthd='LBJ with LU solver';
     case 2
-        mthd='GMRES on stacked system';
+        mthd=sprintf('GMRES on stacked system with ''%s'' preconditioner', options_.simul.preconditioner);
     case 3
-        mthd='BiCGStab on stacked system';
+        mthd=sprintf('BiCGStab on stacked system with ''%s'' preconditioner', options_.simul.preconditioner);
     case 4
         mthd='Sparse LU solver with optimal path length on stacked system';
     case 7
diff --git a/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m b/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
index e3934e2579cefe5092d4b7f6f0bf59559cf56302..67a36b559ad0e67f55874348199c8d94b7016fa2 100644
--- a/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
+++ b/matlab/perfect-foresight-models/solve_two_boundaries_stacked.m
@@ -25,7 +25,7 @@ function [y, T, success, max_res, iter] = solve_two_boundaries_stacked(fh, y, x,
 % ALGORITHM
 %   Newton with LU or GMRES or BiCGStab
 
-% Copyright © 1996-2024 Dynare Team
+% Copyright © 1996-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -57,11 +57,8 @@ verbose = options_.verbosity;
 cvg=false;
 iter=0;
 correcting_factor=0.01;
-ilu_setup.droptol=1e-10;
-ilu_setup.type = 'ilutp';
 max_resa=1e100;
 lambda = 1; % Length of Newton step (unused for stack_solve_algo=4)
-reduced = 0;
 while ~(cvg || iter > options_.simul.maxit)
     r = NaN(Blck_size, periods);
     g1a = spalloc(Blck_size*periods, Blck_size*periods, M_.block_structure.block(Block_Num).NNZDerivatives*periods);
@@ -81,7 +78,6 @@ while ~(cvg || iter > options_.simul.maxit)
             g1a((it_-y_kmin-1)*Blck_size+(1:Blck_size), (it_-y_kmin-2)*Blck_size+(1:3*Blck_size)) = g1;
         end
     end
-    preconditioner = 2;
     ya = reshape(y(y_index, y_kmin+(1:periods)), 1, periods*Blck_size)';
     ra = reshape(r, periods*Blck_size, 1);
     b=-ra+g1a*ya;
@@ -98,7 +94,7 @@ while ~(cvg || iter > options_.simul.maxit)
     end
     if ~cvg
         if iter>0
-            if ~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1)
+            if ~isreal(max_res) || isnan(max_res) %|| (max_resa<max_res && iter>1)
                 if verbose && ~isreal(max_res)
                     disp(['Variable ' M_.endo_names{max_indx} ' (' int2str(max_indx) ') returns an undefined value']);
                 end
@@ -128,7 +124,6 @@ while ~(cvg || iter > options_.simul.maxit)
                     end
                 elseif lambda>1e-8 && stack_solve_algo ~= 4
                     lambda=lambda/2;
-                    reduced = 1;
                     if verbose
                         disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
                     end
@@ -155,102 +150,31 @@ while ~(cvg || iter > options_.simul.maxit)
         g1aa=g1a;
         ba=b;
         max_resa=max_res;
-        if stack_solve_algo==0
+        if stack_solve_algo==0 || (ismember(stack_solve_algo, [2, 3]) && strcmp(options_.simul.preconditioner, 'iterstack') && size(g1a, 1) < options_.simul.iterstack.maxlu) % Fallback to LU if block too small for iterstack
             dx = g1a\b- ya;
             ya = ya + lambda*dx;
             y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
-        elseif stack_solve_algo==2
-            flag1=1;
-            while flag1>0
-                if preconditioner==2
-                    [L1, U1]=ilu(g1a,ilu_setup);
-                elseif preconditioner==3
-                    Size = Blck_size;
-                    gss1 =  g1a(Size + 1: 2*Size,Size + 1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
-                    [L1, U1]=lu(gss1);
-                    L(1:Size,1:Size) = L1;
-                    U(1:Size,1:Size) = U1;
-                    gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
-                    [L2, U2]=lu(gss2);
-                    L(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), L2);
-                    U(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), U2);
-                    gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size);
-                    [L3, U3]=lu(gss2);
-                    L((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = L3;
-                    U((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = U3;
-                    L1 = L;
-                    U1 = U;
-                elseif preconditioner==4
-                    Size = Blck_size;
-                    gss1 =  g1a(1: 3*Size, 1: 3*Size);
-                    [L, U] = lu(gss1);
-                    L1 = kron(eye(ceil(periods/3)),L);
-                    U1 = kron(eye(ceil(periods/3)),U);
-                    L1 = L1(1:periods * Size, 1:periods * Size);
-                    U1 = U1(1:periods * Size, 1:periods * Size);
-                end
-                [za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);
-                if (flag1>0 || reduced)
-                    if verbose
-                        if flag1==1
-                            disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
-                        elseif flag1==2
-                            disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
-                        elseif flag1==3
-                            disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
-                        end
-                    end
-                    ilu_setup.droptol = ilu_setup.droptol/10;
-                    reduced = 0;
-                else
-                    dx = za - ya;
-                    ya = ya + lambda*dx;
-                    y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
-                end
+        elseif ismember(stack_solve_algo, [2, 3])
+            if strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0
+                [L, U, P, Q] = lu(g1a);
+                zc = L\(P*b);
+                zb = U\zc;
+            elseif strcmp(options_.simul.preconditioner, 'iterstack')
+                [L, U, P, Q] = iterstack_preconditioner(g1a, options_);
+            elseif strcmp(options_.simul.preconditioner, 'ilu')
+                [L, U, P] = ilu(g1a, options_.simul.ilu);
+                Q = speye(size(g1a));
             end
-        elseif stack_solve_algo==3
-            flag1=1;
-            while flag1>0
-                if preconditioner==2
-                    [L1, U1]=ilu(g1a,ilu_setup);
-                    [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
-                elseif preconditioner==3
-                    Size = Blck_size;
-                    gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
-                    [L1, U1]=lu(gss0);
-                    P1 = eye(size(gss0));
-                    Q1 = eye(size(gss0));
-                    L = kron(eye(periods),L1);
-                    U = kron(eye(periods),U1);
-                    P = kron(eye(periods),P1);
-                    Q = kron(eye(periods),Q1);
-                    [za,flag1] = bicgstab1(g1a,b,1e-7,Blck_size*periods,L,U, P, Q);
+            if ~(strcmp(options_.simul.preconditioner, 'umfiter') && iter == 0)
+                if stack_solve_algo == 2
+                    zb = gmres(P*g1a*Q, P*b, options_.simul.gmres_restart, options_.simul.itertol, options_.simul.itermaxit, L, U);
                 else
-                    Size = Blck_size;
-                    gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
-                    [L1, U1]=lu(gss0);
-                    L1 = kron(eye(periods),L1);
-                    U1 = kron(eye(periods),U1);
-                    [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
-                end
-                if flag1>0 || reduced
-                    if verbose
-                        if flag1==1
-                            disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
-                        elseif flag1==2
-                            disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
-                        elseif flag1==3
-                            disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
-                        end
-                    end
-                    ilu_setup.droptol = ilu_setup.droptol/10;
-                    reduced = 0;
-                else
-                    dx = za - ya;
-                    ya = ya + lambda*dx;
-                    y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
+                    zb = bicgstab(P*g1a*Q, P*b, options_.simul.itertol, options_.simul.itermaxit, L, U);
                 end
             end
+            dx = Q*zb - ya;
+            ya = ya + lambda*dx;
+            y(y_index, y_kmin+(1:periods))=reshape(ya',length(y_index),periods);
         elseif stack_solve_algo==4
             stpmx = 100 ;
             stpmax = stpmx*max([sqrt(ya'*ya);size(y_index,2)]);