diff --git a/matlab/block_mfs_steadystate.m b/matlab/block_mfs_steadystate.m
index 7496961bd02697b7a196278b8bc6525202086ae0..ceccf3555f0ad8c50f43a3b5767316abeec085b2 100644
--- a/matlab/block_mfs_steadystate.m
+++ b/matlab/block_mfs_steadystate.m
@@ -21,5 +21,5 @@ function [r, g1] = block_mfs_steadystate(y, b, y_all, exo, params, T, M)
 
 y_all(M.block_structure_stat.block(b).variable) = y;
 
-eval(['[r,g1] = ' M.fname '.static(b, y_all, exo, params, T);']);
+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 5f0835879b56c8036ffc46ace22e6784e29acc8a..f1316b28119afeff33ae0b4638285f9037401825 100644
--- a/matlab/dynare_solve_block_or_bytecode.m
+++ b/matlab/dynare_solve_block_or_bytecode.m
@@ -37,18 +37,18 @@ if options.block && ~options.bytecode
                 ss(M.block_structure_stat.block(b).variable) = y;
             else
                 n = length(M.block_structure_stat.block(b).variable);
-                [ss, check] = solve_one_boundary([M.fname '.block.static_' int2str(b)], ss, exo, ...
-                                                 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);
+                [ss, T, check] = solve_one_boundary([M.fname '.block.static_' int2str(b)], ss, exo, ...
+                                                    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
                     info = 1;
                     return
                 end
             end
         end
-        % The following updates the temporary terms vector
-        [r, g1, x, ~, T] = feval([M.fname '.static'], b, ss, exo, params, T);
+        % In particular, the following updates the temporary terms vector (for the dynare_solve case)
+        [r, x, T, g1] = feval([M.fname '.static'], b, ss, exo, params, T);
     end
 elseif options.bytecode
     if options.solve_algo > 4
diff --git a/matlab/evaluate_static_model.m b/matlab/evaluate_static_model.m
index f418d122d2f456c84b287d78c0bd2c698d181166..a4aa46ddeaf8a68e56c6526c8b4667e100259824 100644
--- a/matlab/evaluate_static_model.m
+++ b/matlab/evaluate_static_model.m
@@ -52,12 +52,12 @@ else
             % have zero residuals by construction
             if M.block_structure_stat.block(b).Simulation_Type ~= 1 && ...
                     M.block_structure_stat.block(b).Simulation_Type ~= 2
-                [r, ~, ~, ~, T] = feval(fh_static,b,ys,exo_ss,params,T);
+                [r, ~, T] = feval(fh_static,b,ys,exo_ss,params,T);
                 residuals(mfsb) = r;
             else
                 %need to evaluate the recursive blocks to compute the
                 %temporary terms
-                [~, ~, ~, ~, T] = feval(fh_static,b,ys,exo_ss,params,T);
+                [~, ~, T] = feval(fh_static,b,ys,exo_ss,params,T);
             end
         end
     else
diff --git a/matlab/lnsrch1_wrapper_one_boundary.m b/matlab/lnsrch1_wrapper_one_boundary.m
index 8ac9c14d0cb8c0d16cf2234d3fd908cfbc273dd3..f10f9c80d089883405ac67f6fa49e47d687c251a 100644
--- a/matlab/lnsrch1_wrapper_one_boundary.m
+++ b/matlab/lnsrch1_wrapper_one_boundary.m
@@ -1,4 +1,4 @@
-function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, steady_state, it_)
+function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, steady_state, T, it_)
 % wrapper for solve_one_boundary m-file when it is used with a dynamic
 % model
 %
@@ -11,6 +11,8 @@ function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, stea
 %   y                   [matrix]        All the endogenous variables of the model
 %   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
 % OUTPUTS
 %   r                   [vector]        The residuals of the current block
 %
@@ -21,7 +23,7 @@ function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, stea
 %   none.
 %
 
-% Copyright (C) 2009-2017 Dynare Team
+% Copyright (C) 2009-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -40,4 +42,4 @@ function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, stea
 
 %reshape the input arguments of the dynamic function
 y(it_, :) = ya;
-[r, y, g1, g2, g3]=feval(fname, y, x, params, steady_state, it_, 0);
+[r, y, T, g1, g2, g3]=feval(fname, y, x, params, steady_state, T, it_, false);
diff --git a/matlab/lnsrch1_wrapper_two_boundaries.m b/matlab/lnsrch1_wrapper_two_boundaries.m
index 05c396c77fcd32f93d4437cd4fe3c8f8509fe9bd..127bf15e0797eabc9ec269894a8d9aac946435cc 100644
--- a/matlab/lnsrch1_wrapper_two_boundaries.m
+++ b/matlab/lnsrch1_wrapper_two_boundaries.m
@@ -1,5 +1,5 @@
 function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
-                                             params, steady_state, periods, y_kmin, y_size,Periods)
+                                             params, steady_state, T, periods, y_kmin, y_size,Periods)
 % wrapper for solve_one_boundary m-file when it is used with a dynamic
 % model
 %
@@ -12,6 +12,8 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
 %   y                   [matrix]        All the endogenous variables of the model
 %   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
 %   periods             [int]           The number of periods
 %   y_kmin              [int]           The maximum number of lag on en endogenous variables
 %   y_size              [int]           The number of endogenous variables
@@ -26,7 +28,7 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
 %   none.
 %
 
-% Copyright (C) 2009-2017 Dynare Team
+% Copyright (C) 2009-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -45,5 +47,5 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
 
 %reshape the input arguments of the dynamic function
 y(y_kmin+1:y_kmin+periods, y_index) = reshape(ya',length(y_index),periods)';
-[r, y, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, periods, 0, y_kmin, y_size, Periods);
+[r, y, T, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, y_size, Periods);
 ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*y_size, 1);
diff --git a/matlab/resid.m b/matlab/resid.m
index 807dc494f799f5aaf728b35a40cc35493fb72ea2..46d876053c9e803476933778e6ce86861e64dbac 100644
--- a/matlab/resid.m
+++ b/matlab/resid.m
@@ -71,7 +71,7 @@ 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, T] = feval([M_.fname '.static'],...
+        [r, yy, T, g, var_indx] = feval([M_.fname '.static'],...
                                         i,...
                                         oo_.steady_state,...
                                         [oo_.exo_steady_state; ...
diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m
index 39a073c225d9d01e0aa4003cec4153b0cf100a4e..b0cf6b5be49d15bdf0aea75b0e4cc8807f0233bd 100644
--- a/matlab/solve_one_boundary.m
+++ b/matlab/solve_one_boundary.m
@@ -1,5 +1,5 @@
-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)
+function [y, T, 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
 %
@@ -42,6 +42,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, T, ..
 %
 % OUTPUTS
 %   y                  [matrix]         All endogenous variables of the model
+%   T                  [matrix]        Temporary terms
 %   info               [integer]        >=0 no error
 %                                       <0 error
 %
@@ -93,9 +94,9 @@ for it_=start:incr:finish
     g1=spalloc( Blck_size, Blck_size, nze);
     while ~(cvg==1 || iter>maxit_)
         if is_dynamic
-            [r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, it_, 0);
+            [r, y, T, g1, g2, g3] = feval(fname, y, x, params, steady_state, T, it_, false);
         else
-            [r, y, g1] = feval(fname, y, x, params, T);
+            [r, y, T, g1] = feval(fname, y, x, params, T);
         end
         if ~isreal(r)
             max_res=(-(max(max(abs(r))))^2)^0.5;
@@ -204,7 +205,9 @@ for it_=start:incr:finish
                 p = -g1\r ;
                 [ya,f,r,check]=lnsrch1(y(it_,:)',f,g,p,stpmax, ...
                                        'lnsrch1_wrapper_one_boundary',nn, ...
-                                       y_index_eq, options.solve_tolx, y_index_eq, fname, y, x, params, steady_state, it_);
+                                       y_index_eq, options.solve_tolx, y_index_eq, fname, y, x, params, steady_state, T, it_);
+                %% Recompute temporary terms, since they are not given as output of lnsrch1
+                [~, ~, T] = feval(fname, y, x, params, steady_state, T, it_, false);
                 dx = ya' - y(it_, :);
                 y(it_,:) = ya';
             elseif (is_dynamic && (stack_solve_algo==1 || stack_solve_algo==0)) || (~is_dynamic && options.solve_algo==6)
@@ -261,10 +264,10 @@ for it_=start:incr:finish
                         y(y_index_eq) = phat;
                     end
                     if is_dynamic
-                        [r, y, g1, g2, g3] = feval(fname, y, x, params, ...
-                                                   steady_state, it_, 0);
+                        [r, y, T, g1, g2, g3] = feval(fname, y, x, params, ...
+                                                   steady_state, T, it_, false);
                     else
-                        [r, y, g1] = feval(fname, y, x, params, T);
+                        [r, y, T, 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/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m
index 4376ab82c0f20b55bee01f6551098ce063a65b1a..9eda54694346d5d2800aea7fc672a6bcc53aa2d8 100644
--- a/matlab/solve_two_boundaries.m
+++ b/matlab/solve_two_boundaries.m
@@ -1,4 +1,4 @@
-function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_index, nze, periods, y_kmin_l, y_kmax_l, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo,options,M, oo)
+function [y, T, oo]= solve_two_boundaries(fname, y, x, params, steady_state, T, y_index, nze, periods, y_kmin_l, y_kmax_l, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo,options,M, oo)
 % Computes the deterministic simulation of a block of equation containing
 % both lead and lag variables using relaxation methods
 %
@@ -9,6 +9,7 @@ function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_inde
 %   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             [vector of int] The index of the endogenous variables of
 %                                       the block
 %   nze                 [integer]       number of non-zero elements in the
@@ -36,6 +37,7 @@ function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_inde
 %
 % OUTPUTS
 %   y                   [matrix]        All endogenous variables of the model
+%   T                   [matrix]        Temporary terms
 %   oo                  [structure]     Results
 %
 % ALGORITHM
@@ -82,7 +84,7 @@ Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods);
 g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
 reduced = 0;
 while ~(cvg==1 || iter>maxit_)
-    [r, y, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, periods, 0, y_kmin, Blck_size,options.periods);
+    [r, y, T, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, Blck_size,options.periods);
     preconditioner = 2;
     g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
     term1 = g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)';
@@ -306,7 +308,7 @@ while ~(cvg==1 || iter>maxit_)
             g = (ra'*g1a)';
             f = 0.5*ra'*ra;
             p = -g1a\ra;
-            [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, options.solve_tolx, fname, y, y_index,x, params, steady_state, periods, y_kmin, Blck_size,options.periods);
+            [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, options.solve_tolx, fname, y, y_index,x, params, steady_state, T, periods, y_kmin, Blck_size,options.periods);
             dx = ya - yn;
             y(1+y_kmin:periods+y_kmin,y_index)=reshape(yn',length(y_index),periods)';
         end
diff --git a/preprocessor b/preprocessor
index 04b7d4386dd8a6dbd74c032e92d09e6ca0758fd6..faa6666abef430f777fd43d175de668611d67260 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 04b7d4386dd8a6dbd74c032e92d09e6ca0758fd6
+Subproject commit faa6666abef430f777fd43d175de668611d67260