diff --git a/matlab/block_mfs_steadystate.m b/matlab/block_mfs_steadystate.m
index a0793e674d99e79e15125a8a530f85d5ef3b80b5..7496961bd02697b7a196278b8bc6525202086ae0 100644
--- a/matlab/block_mfs_steadystate.m
+++ b/matlab/block_mfs_steadystate.m
@@ -1,8 +1,8 @@
-function [r, g1] = block_mfs_steadystate(y, b, y_all, exo, params, M)
+function [r, g1] = block_mfs_steadystate(y, b, y_all, exo, params, T, M)
 % Wrapper around the *_static.m file, for use with dynare_solve,
 % when block_mfs option is given to steady.
 
-% Copyright (C) 2009-2012 Dynare Team
+% Copyright (C) 2009-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -21,5 +21,5 @@ function [r, g1] = block_mfs_steadystate(y, b, y_all, exo, params, M)
 
 y_all(M.block_structure_stat.block(b).variable) = y;
 
-eval(['[r,g1] = ' M.fname '.static(b, y_all, exo, params);']);
+eval(['[r,g1] = ' M.fname '.static(b, y_all, exo, params, T);']);
 g1 = full(g1);
diff --git a/matlab/dynare_solve_block_or_bytecode.m b/matlab/dynare_solve_block_or_bytecode.m
index 69cc8d895787600a7ab54edda3b5e76d40a62046..5f0835879b56c8036ffc46ace22e6784e29acc8a 100644
--- a/matlab/dynare_solve_block_or_bytecode.m
+++ b/matlab/dynare_solve_block_or_bytecode.m
@@ -19,6 +19,7 @@ function [x,info] = dynare_solve_block_or_bytecode(y, exo, params, options, M)
 info = 0;
 x = y;
 if options.block && ~options.bytecode
+    T = NaN(M.block_structure_stat.tmp_nbr, 1);
     for b = 1:length(M.block_structure_stat.block)
         ss = x;
         if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
@@ -26,7 +27,7 @@ if options.block && ~options.bytecode
             if options.solve_algo <= 4 || options.solve_algo >= 9
                 [y, check] = dynare_solve('block_mfs_steadystate', ...
                                           ss(M.block_structure_stat.block(b).variable), ...
-                                          options, b, ss, exo, params, M);
+                                          options, b, ss, exo, params, T, M);
                 if check ~= 0
                     %                    error(['STEADY: convergence
                     %                    problems in block ' int2str(b)])
@@ -37,7 +38,7 @@ if options.block && ~options.bytecode
             else
                 n = length(M.block_structure_stat.block(b).variable);
                 [ss, check] = solve_one_boundary([M.fname '.block.static_' int2str(b)], ss, exo, ...
-                                                 params, [], M.block_structure_stat.block(b).variable, n, 1, false, b, 0, options.simul.maxit, ...
+                                                 params, [], T, M.block_structure_stat.block(b).variable, n, 1, false, b, 0, options.simul.maxit, ...
                                                  options.solve_tolf, ...
                                                  options.slowc, 0, options.solve_algo, true, false, false, M, options);
                 if check
@@ -46,8 +47,8 @@ if options.block && ~options.bytecode
                 end
             end
         end
-        [r, g1, x] = feval([M.fname '.static'], b, ss, ...
-                           exo, params);
+        % The following updates the temporary terms vector
+        [r, g1, x, ~, T] = feval([M.fname '.static'], b, ss, exo, params, T);
     end
 elseif options.bytecode
     if options.solve_algo > 4
diff --git a/matlab/resid.m b/matlab/resid.m
index 1f3743a072e7c590931eb1e9f0986414b672738a..807dc494f799f5aaf728b35a40cc35493fb72ea2 100644
--- a/matlab/resid.m
+++ b/matlab/resid.m
@@ -69,12 +69,13 @@ end
 % Compute the residuals
 if options_.block && ~options_.bytecode
     z = zeros(M_.endo_nbr,1);
+    T = NaN(M_.block_structure_stat.tmp_nbr, 1);
     for i = 1:length(M_.block_structure_stat.block)
-        [r, g, yy, var_indx] = feval([M_.fname '.static'],...
-                                     i,...
-                                     oo_.steady_state,...
-                                     [oo_.exo_steady_state; ...
-                            oo_.exo_det_steady_state], M_.params);
+        [r, g, yy, var_indx, T] = feval([M_.fname '.static'],...
+                                        i,...
+                                        oo_.steady_state,...
+                                        [oo_.exo_steady_state; ...
+                                         oo_.exo_det_steady_state], M_.params, T);
         idx = M_.block_structure_stat.block(i).equation;
         z(idx) = r;
     end
diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m
index 503131c46101c6e0cae80bdca4c32d05863d484c..39a073c225d9d01e0aa4003cec4153b0cf100a4e 100644
--- a/matlab/solve_one_boundary.m
+++ b/matlab/solve_one_boundary.m
@@ -1,4 +1,4 @@
-function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ...
+function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, T, ...
                                         y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo, is_forward, is_dynamic, verbose, M, options, oo)
 % Computes the deterministic simulation of a block of equation containing
 % lead or lag variables
@@ -10,6 +10,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ...
 %   x                   [matrix]        All the exogenous variables of the model
 %   params              [vector]        All the parameters of the model
 %   steady_state        [vector]        steady state of the model
+%   T                   [matrix]        Temporary terms
 %   y_index_eq          [vector of int] The index of the endogenous variables of
 %                                       the block
 %   nze                 [integer]       number of non-zero elements in the
@@ -94,7 +95,7 @@ for it_=start:incr:finish
         if is_dynamic
             [r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, it_, 0);
         else
-            [r, y, g1] = feval(fname, y, x, params);
+            [r, y, g1] = feval(fname, y, x, params, T);
         end
         if ~isreal(r)
             max_res=(-(max(max(abs(r))))^2)^0.5;
@@ -263,7 +264,7 @@ for it_=start:incr:finish
                         [r, y, g1, g2, g3] = feval(fname, y, x, params, ...
                                                    steady_state, it_, 0);
                     else
-                        [r, y, g1] = feval(fname, y, x, params);
+                        [r, y, g1] = feval(fname, y, x, params, T);
                     end
                     if max(abs(r))>=options.solve_tolf
                         [dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
diff --git a/preprocessor b/preprocessor
index a6d9ba6e558c9f149ee81cc0e00837623563cb50..04b7d4386dd8a6dbd74c032e92d09e6ca0758fd6 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit a6d9ba6e558c9f149ee81cc0e00837623563cb50
+Subproject commit 04b7d4386dd8a6dbd74c032e92d09e6ca0758fd6