From c6bac390d45bb8f3efa94123c977766c940700b0 Mon Sep 17 00:00:00 2001
From: Houtan Bastani <houtan@dynare.org>
Date: Fri, 21 Apr 2017 18:15:32 +0200
Subject: [PATCH] add exogenous steady state to dynamic function call. #825

---
 matlab/discretionary_policy_1.m                |  5 +++--
 matlab/dyn_risky_steadystate_solver.m          | 17 +++++++++++------
 matlab/evaluate_max_dynamic_residual.m         |  6 +++---
 matlab/evaluate_steady_state.m                 |  7 ++++---
 matlab/getH.m                                  | 10 +++++-----
 matlab/identification_analysis.m               |  4 ++--
 matlab/model_diagnostics.m                     |  6 +++---
 matlab/partial_information/dr1_PI.m            |  5 +++--
 .../det_cond_forecast.m                        |  8 +++++---
 .../perfect_foresight_mcp_problem.m            | 11 ++++++-----
 .../perfect_foresight_problem.m                | 11 ++++++-----
 .../perfect_foresight_solver_core.m            | 18 +++++++++---------
 .../private/simulation_core.m                  | 16 ++++++++--------
 matlab/perfect-foresight-models/sim1.m         |  6 +++---
 matlab/perfect-foresight-models/sim1_lbj.m     |  8 ++++----
 matlab/perfect-foresight-models/sim1_linear.m  |  6 +++---
 .../sim1_purely_backward.m                     |  6 +++---
 .../sim1_purely_forward.m                      |  6 +++---
 .../solve_stacked_linear_problem.m             |  4 ++--
 .../solve_stacked_problem.m                    |  6 +++---
 matlab/simul_backward_linear_model.m           |  4 ++--
 matlab/simul_backward_nonlinear_model.m        |  4 ++--
 matlab/stochastic_solvers.m                    |  7 ++++---
 matlab/thet2tau.m                              |  4 ++--
 tests/lmmcp/sw_newton.mod                      |  4 ++--
 25 files changed, 101 insertions(+), 88 deletions(-)

diff --git a/matlab/discretionary_policy_1.m b/matlab/discretionary_policy_1.m
index 931d3b596f..7a987f7b8d 100644
--- a/matlab/discretionary_policy_1.m
+++ b/matlab/discretionary_policy_1.m
@@ -1,6 +1,6 @@
 function [dr,ys,info]=discretionary_policy_1(oo_,Instruments)
 
-% Copyright (C) 2007-2016 Dynare Team
+% Copyright (C) 2007-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -86,7 +86,8 @@ if exo_nbr == 0
 end
 
 [junk,jacobia_] = feval([M_.fname '_dynamic'],z, [zeros(size(oo_.exo_simul)) ...
-                    oo_.exo_det_simul], M_.params, zeros(endo_nbr,1), it_);
+                    oo_.exo_det_simul], M_.params, zeros(endo_nbr,1), ...
+                        zeros(exo_nbr,1), it_);
 if any(junk~=0)
     error(['discretionary_policy: the model must be written in deviation ' ...
            'form and not have constant terms'])
diff --git a/matlab/dyn_risky_steadystate_solver.m b/matlab/dyn_risky_steadystate_solver.m
index a7e09b1bd3..7864b86f97 100644
--- a/matlab/dyn_risky_steadystate_solver.m
+++ b/matlab/dyn_risky_steadystate_solver.m
@@ -62,7 +62,7 @@ function [dr,info] = dyn_risky_steadystate_solver(ys0,M, ...
     %! @end deftypefn
     %@eod:
 
-    % Copyright (C) 2001-2012 Dynare Team
+    % Copyright (C) 2001-2017 Dynare Team
     %
     % This file is part of Dynare.
     %
@@ -162,7 +162,8 @@ function [resid,dr] = risky_residuals(ys,pm,M,dr,options,oo)
     z = z(iyr0) ;
     [resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
                            [oo.exo_simul ...
-                        oo.exo_det_simul], M.params, dr.ys, 2);
+                        oo.exo_det_simul], M.params, dr.ys, oo.exo_steady_state, ...
+                           2);
     if ~isreal(d1) || ~isreal(d2)
         pause
     end
@@ -222,7 +223,8 @@ function [resid,dr] = risky_residuals_ds(x,pm,M,dr,options,oo)
     z = z(iyr0) ;
     [resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
                            [oo.exo_simul ...
-                        oo.exo_det_simul], M.params, dr.ys, 2);
+                        oo.exo_det_simul], M.params, dr.ys, oo.exo_steady_state, ...
+                           2);
     if ~isreal(d1) || ~isreal(d2)
         pause
     end
@@ -255,7 +257,8 @@ function dr_np = first_step_ds(x,pm,M,dr,options,oo)
     z = z(iyr0) ;
     [resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
                            [oo.exo_simul ...
-                        oo.exo_det_simul], M.params, dr.ys, 2);
+                        oo.exo_det_simul], M.params, dr.ys, oo.exo_steady_state, ...
+                           2);
     if ~isreal(d1) || ~isreal(d2)
         pause
     end
@@ -296,7 +299,8 @@ function [resid,dr] = risky_residuals_k_order(ys,pm,M,dr,options,oo)
     z = z(iyr0) ;
     [resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
                               [oo.exo_simul ...
-                        oo.exo_det_simul], M.params, dr.ys, 2);
+                        oo.exo_det_simul], M.params, dr.ys, oo.exo_steady_state, ...
+                           2);
     
     if isfield(options,'portfolio') && options.portfolio == 1
         eq_np = pm.eq_np;
@@ -331,7 +335,8 @@ function [resid,dr] = risky_residuals_k_order(ys,pm,M,dr,options,oo)
     if nargout > 1
         [resid1,d1,d2,d3] = feval([M.fname '_dynamic'],z,...
                                   [oo.exo_simul ...
-                            oo.exo_det_simul], M.params, dr.ys, 2);
+                            oo.exo_det_simul], M.params, dr.ys, oo.exo_steady_state, ...
+                                  2);
 
         
         [a,b,c] = find(d2(eq_np,pm.i_d2_np));
diff --git a/matlab/evaluate_max_dynamic_residual.m b/matlab/evaluate_max_dynamic_residual.m
index 41aa4ec249..aea6167bba 100644
--- a/matlab/evaluate_max_dynamic_residual.m
+++ b/matlab/evaluate_max_dynamic_residual.m
@@ -1,6 +1,6 @@
-function err = evaluate_max_dynamic_residual(model_dynamic, Y, exogenous_variables, params, steady_state, periods, ny, max_lag, lead_lag_incidence)
+function err = evaluate_max_dynamic_residual(model_dynamic, Y, exogenous_variables, params, steady_state, steady_state_x, periods, ny, max_lag, lead_lag_incidence)
 
-% Copyright (C) 2013 Dynare Team
+% Copyright (C) 2013-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -23,7 +23,7 @@ i_cols = find(lead_lag_incidence');
 err = 0;
 
 for it = (max_lag+1):(max_lag+periods)
-    d = model_dynamic(Y(i_cols), exogenous_variables, params, steady_state, it);
+    d = model_dynamic(Y(i_cols), exogenous_variables, params, steady_state, steady_state_x, it);
     i_rows = i_rows + ny;
     i_cols = i_cols + ny;
     r = max(abs(d));
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 053d62255e..c2ef05920d 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,M,options,oo,steadysta
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2001-2016 Dynare Team
+% Copyright (C) 2001-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -313,12 +313,13 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
             [chck, r, junk]= bytecode('dynamic','evaluate', z, zx, M.params, ys, 1);
             mexErrCheck('bytecode', chck);
         elseif options.block
-            [r, oo.dr] = feval([M.fname '_dynamic'], z', zx, M.params, ys, M.maximum_lag+1, oo.dr);
+            [r, oo.dr] = feval([M.fname '_dynamic'], z', zx, M.params, ys, ...
+                               M.maximum_lag+1, oo.dr);
         else
             iyv = M.lead_lag_incidence';
             iyr0 = find(iyv(:));
             xys = z(iyr0);
-            r = feval([M.fname '_dynamic'], z(iyr0), zx, M.params, ys, M.maximum_lag + 1);
+            r = feval([M.fname '_dynamic'], z(iyr0), zx, M.params, ys, oo.exo_steady_state, M.maximum_lag + 1);
         end
 
         % Fail if residual greater than tolerance
diff --git a/matlab/getH.m b/matlab/getH.m
index 0db11d8e33..d114f152ec 100644
--- a/matlab/getH.m
+++ b/matlab/getH.m
@@ -32,7 +32,7 @@ function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, estim_params_,M_,oo
 %   H2s:            Hessian of steady state with respect to estimated
 %                   structural parameters only (indx)
 
-% Copyright (C) 2010-2016 Dynare Team
+% Copyright (C) 2010-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -125,7 +125,7 @@ end
 if kronflag==-2,
     if nargout>5,
         [residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
-            M_.params, oo_.dr.ys, 1);
+            M_.params, oo_.dr.ys, oo_.exo_steady_state, 1);
         g22 = hessian_sparse('thet2tau',[M_.params(indx)],options_.gstep,estim_params_,M_, oo_, indx,[],-1);
         H2ss=full(g22(1:M_.endo_nbr,:));
         H2ss = reshape(H2ss,[M_.endo_nbr param_nbr param_nbr]);
@@ -145,7 +145,7 @@ if kronflag==-2,
         clear gx22;
     else
         [residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
-            M_.params, oo_.dr.ys, 1);        
+            M_.params, oo_.dr.ys, oo_.exo_steady_state, 1);
     end
     gp = fjaco('thet2tau',[M_.params(indx)],estim_params_,M_, oo_, indx,[],-1);
     Hss=gp(1:M_.endo_nbr,:);
@@ -166,7 +166,7 @@ dyssdtheta = -gg1\df;
 if nargout>5,
     [residual, gg1, gg2] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
     [residual, g1, g2, g3] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
-        M_.params, oo_.dr.ys, 1);
+        M_.params, oo_.dr.ys, oo_.exo_steady_state, 1);
     [nr, nc]=size(gg2);
 
     [df, gpx, d2f] = feval([M_.fname,'_static_params_derivs'],oo_.dr.ys, oo_.exo_steady_state', ...
@@ -218,7 +218,7 @@ else
     [df, gp] = feval([M_.fname,'_params_derivs'],yy0, repmat(oo_.exo_steady_state',[M_.maximum_exo_lag+M_.maximum_exo_lead+1,1]), ...
         M_.params, oo_.dr.ys, 1, dyssdtheta,d2yssdtheta);
     [residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, repmat(oo_.exo_steady_state',[M_.maximum_exo_lag+M_.maximum_exo_lead+1,1]), ...
-        M_.params, oo_.dr.ys, 1);
+        M_.params, oo_.dr.ys, oo_.exo_steady_state, 1);
 end
 
 [nr, nc]=size(g2);
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index 88b60ca4d4..6c2c085433 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -28,7 +28,7 @@ function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info, opt
 % SPECIAL REQUIREMENTS
 %    None
 
-% Copyright (C) 2008-2016 Dynare Team
+% Copyright (C) 2008-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -79,7 +79,7 @@ if info(1)==0,
     yy0=oo_.dr.ys(I);
     [residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, ...
         repmat(oo_.exo_steady_state',[M_.maximum_exo_lag+M_.maximum_exo_lead+1]), M_.params, ...
-        oo_.dr.ys, 1);
+        oo_.dr.ys, oo_.exo_steady_state, 1);
     vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];
 
     [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m
index 32a01725d5..b96ec7d48d 100644
--- a/matlab/model_diagnostics.m
+++ b/matlab/model_diagnostics.m
@@ -16,7 +16,7 @@ function model_diagnostics(M,options,oo)
 %   none.
 %  
 
-% Copyright (C) 1996-2013 Dynare Team
+% Copyright (C) 1996-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -209,7 +209,7 @@ if ~options.block
             jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd];
         else
             [junk,jacobia_] = feval([M.fname '_dynamic'],z(iyr0),exo_simul, ...
-                M.params, dr.ys, it_);
+                M.params, dr.ys, oo.exo_steady_state, it_);
         end;
     elseif options.order >= 2
         if (options.bytecode)
@@ -219,7 +219,7 @@ if ~options.block
         else
             [junk,jacobia_,hessian1] = feval([M.fname '_dynamic'],z(iyr0),...
                 exo_simul, ...
-                M.params, dr.ys, it_);
+                M.params, dr.ys, oo.exo_steady_state, it_);
         end;
         if options.use_dll
             % In USE_DLL mode, the hessian is in the 3-column sparse representation
diff --git a/matlab/partial_information/dr1_PI.m b/matlab/partial_information/dr1_PI.m
index 0475d4befc..76c9405921 100644
--- a/matlab/partial_information/dr1_PI.m
+++ b/matlab/partial_information/dr1_PI.m
@@ -40,7 +40,7 @@ function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
 %   none.
 %  
 
-% Copyright (C) 1996-2012 Dynare Team
+% Copyright (C) 1996-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -131,7 +131,8 @@ else
     end
     z = z(iyr0) ;
     [junk,jacobia_] = feval([M_.fname '_dynamic'],z,[oo_.exo_simul ...
-                        oo_.exo_det_simul], M_.params, dr.ys, it_);
+                        oo_.exo_det_simul], M_.params, dr.ys, oo_.exo_steady_state, ...
+                            it_);
 
     if options_.ACES_solver==1 && (length(sim_ruleids)>0 || length(tct_ruleids)>0 )
         if length(sim_ruleids)>0
diff --git a/matlab/perfect-foresight-models/det_cond_forecast.m b/matlab/perfect-foresight-models/det_cond_forecast.m
index e261de89c3..b21d262665 100644
--- a/matlab/perfect-foresight-models/det_cond_forecast.m
+++ b/matlab/perfect-foresight-models/det_cond_forecast.m
@@ -12,7 +12,7 @@ function data_set = det_cond_forecast(varargin)
 %  dataset                [dseries]     Returns a dseries containing the forecasted endgenous variables and shocks
 %
 %
-% Copyright (C) 2013-2016 Dynare Team
+% Copyright (C) 2013-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -475,7 +475,8 @@ if pf && ~surprise
                 if (options_.bytecode)
                     [chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
                 else
-                    [zz, g1b] = feval([M_.fname '_dynamic'], z', zx, M_.params, oo_.steady_state, k);
+                    [zz, g1b] = feval([M_.fname '_dynamic'], z', zx, M_.params, ...
+                                      oo_.steady_state, oo_.exo_steady_state, k);
                     data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end);
                     data1.g1 = g1b(:,1 : end - M_.exo_nbr);
                     chck = 0;
@@ -747,7 +748,8 @@ else
                     if (options_.bytecode)
                         [chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
                     else
-                        [zz, g1b] = feval([M_.fname '_dynamic'], z', zx, M_.params, oo_.steady_state, k);
+                        [zz, g1b] = feval([M_.fname '_dynamic'], z', zx, ...
+                                          M_.params, oo_.steady_state, oo_.exo_steady_state, k);
                         data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end);
                         data1.g1 = g1b(:,1 : end - M_.exo_nbr);
                         chck = 0;
diff --git a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
index d4b4193086..fbf0fa6494 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
@@ -1,10 +1,10 @@
 function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_function, Y0, YT, ...
-                                           exo_simul, params, steady_state, ...
+                                           exo_simul, params, steady_state, exo_steady_state, ...
                                            maximum_lag, T, ny, i_cols, ...
                                            i_cols_J1, i_cols_1, i_cols_T, ...
                                            i_cols_j,nnzJ,eq_index)
 % function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_function, Y0, YT, ...
-%                                            exo_simul, params, steady_state, ...
+%                                            exo_simul, params, steady_state, exo_steady_state, ...
 %                                            maximum_lag, T, ny, i_cols, ...
 %                                            i_cols_J1, i_cols_1, i_cols_T, ...
 %                                            i_cols_j,nnzJ,eq_index)
@@ -20,6 +20,7 @@ function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_functi
 %                                for all simulation periods           
 %   params              [double] nparams*1 array, parameter values
 %   steady_state        [double] endo_nbr*1 vector of steady state values
+%   exo_steady_state    [double] exo_nbr*1 vector of steady state values
 %   maximum_lag         [scalar] maximum lag present in the model
 %   T                   [scalar] number of simulation periods
 %   ny                  [scalar] number of endogenous variables
@@ -45,7 +46,7 @@ function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_functi
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright (C) 1996-2016 Dynare Team
+% Copyright (C) 1996-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -77,11 +78,11 @@ i_cols_J = i_cols;
 for it = 2:(T+1)
     if nargout == 1
         res = dynamic_function(YY(i_cols),exo_simul, params, ...
-            steady_state,it);
+            steady_state, exo_steady_state, it);
         residuals(i_rows) = res(eq_index);
     elseif nargout == 2
         [res,jacobian] = dynamic_function(YY(i_cols),exo_simul, params, ...
-            steady_state,it);
+            steady_state, exo_steady_state, it);
         residuals(i_rows) = res(eq_index);
         if it == 2
             [rows,cols,vals] = find(jacobian(eq_index,i_cols_1));
diff --git a/matlab/perfect-foresight-models/perfect_foresight_problem.m b/matlab/perfect-foresight-models/perfect_foresight_problem.m
index 9561190945..311eac3616 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_problem.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_problem.m
@@ -1,10 +1,10 @@
 function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, Y0, YT, ...
-                                           exo_simul, params, steady_state, ...
+                                           exo_simul, params, steady_state, exo_steady_state, ...
                                            maximum_lag, T, ny, i_cols, ...
                                            i_cols_J1, i_cols_1, i_cols_T, ...
                                            i_cols_j,nnzJ)
 % function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, Y0, YT, ...
-%                                            exo_simul, params, steady_state, ...
+%                                            exo_simul, params, steady_state, exo_steady_state, ...
 %                                            maximum_lag, T, ny, i_cols, ...
 %                                            i_cols_J1, i_cols_1, i_cols_T, ...
 %                                            i_cols_j,nnzJ)
@@ -19,6 +19,7 @@ function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function,
 %                                for all simulation periods           
 %   params              [double] nparams*1 array, parameter values
 %   steady_state        [double] endo_nbr*1 vector of steady state values
+%   exo_steady_state    [double] exo_nbr*1 vector of steady state values
 %   maximum_lag         [scalar] maximum lag present in the model
 %   T                   [scalar] number of simulation periods
 %   ny                  [scalar] number of endogenous variables
@@ -42,7 +43,7 @@ function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function,
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright (C) 1996-2016 Dynare Team
+% Copyright (C) 1996-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -74,10 +75,10 @@ function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function,
     for it = 2:(T+1)
         if nargout == 1
              residuals(i_rows) = dynamic_function(YY(i_cols),exo_simul, params, ...
-                                                         steady_state,it);
+                                                         steady_state, exo_steady_state, it);
         elseif nargout == 2
             [residuals(i_rows),jacobian] = dynamic_function(YY(i_cols),exo_simul, params, ...
-                                                         steady_state,it);
+                                                         steady_state, exo_steady_state, it);
             if it == 2
                 [rows,cols,vals] = find(jacobian(:,i_cols_1));
                 iJacobian{1} = [offset+rows, i_cols_J1(cols), vals];
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index da267b7dae..d3bd13a1de 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -11,7 +11,7 @@ function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_)
 % - oo_                 [struct] contains results
 % - maxerror            [double] contains the maximum absolute error
 
-% Copyright (C) 2015-2016 Dynare Team
+% Copyright (C) 2015-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -79,36 +79,36 @@ else
     else
         if M_.maximum_endo_lead == 0 % Purely backward model
             [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                sim1_purely_backward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                sim1_purely_backward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
         elseif M_.maximum_endo_lag == 0 % Purely forward model
             [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                sim1_purely_forward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                sim1_purely_forward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
         else % General case
             switch options_.stack_solve_algo
               case 0
                 if options_.linear_approximation
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                        sim1_linear(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
+                        sim1_linear(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_, oo_);
                 else
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                        sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                        sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
                 end
               case 6
                 if options_.linear_approximation
                     error('Invalid value of stack_solve_algo option!')
                 end
                 [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                    sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                    sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
               case 7
                 if options_.linear_approximation
                     if isequal(options_.solve_algo, 10)
                         warning('It would be more efficient to set option solve_algo equal to 0!')
                     end
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                        solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
+                        solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_, oo_);
                 else
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                        solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                        solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
                 end
               otherwise
                 error('Invalid value of stack_solve_algo option!')
@@ -135,7 +135,7 @@ if nargout>1
             [chck, residuals, junk]= bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
         else
             residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '_dynamic']), y0, yT, ...
-                                                  oo_.exo_simul,M_.params,oo_.steady_state, ...
+                                                  oo_.exo_simul,M_.params,oo_.steady_state,oo_.exo_steady_state, ...
                                                   M_.maximum_lag,options_.periods,M_.endo_nbr,i_cols, ...
                                                   i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ...
                                                   M_.NNZDerivatives(1));
diff --git a/matlab/perfect-foresight-models/private/simulation_core.m b/matlab/perfect-foresight-models/private/simulation_core.m
index b52ecad172..0afdba517e 100644
--- a/matlab/perfect-foresight-models/private/simulation_core.m
+++ b/matlab/perfect-foresight-models/private/simulation_core.m
@@ -1,7 +1,7 @@
 function [oo_, maxerror] = simulation_core(options_, M_, oo_)
 %function [oo_, maxerror] = simulation_core(options_, M_, oo_)
 
-% Copyright (C) 2015 Dynare Team
+% Copyright (C) 2015-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -64,29 +64,29 @@ else
     else
         if M_.maximum_endo_lead == 0 % Purely backward model
             [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                sim1_purely_backward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                sim1_purely_backward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
         elseif M_.maximum_endo_lag == 0 % Purely forward model
             [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                sim1_purely_forward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                sim1_purely_forward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
         else % General case
             if options_.stack_solve_algo == 0
                 if options_.linear_approximation
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                        sim1_linear(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
+                        sim1_linear(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_, oo_);
                 else
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                        sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                        sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
                 end
             elseif options_.stack_solve_algo == 6
                 [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                    sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                    sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
             elseif options_.stack_solve_algo == 7
                 if options_.linear_approximation
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                        solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
+                        solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_, oo_);
                 else
                     [oo_.endo_simul, oo_.deterministic_simulation] = ...
-                        solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
+                        solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_, oo_);
                 end
             end
         end
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 1a173b2115..e1a6d42812 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -1,4 +1,4 @@
-function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options)
+function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options, oo)
 
 % Performs deterministic simulations with lead or lag on one period. Uses sparse matrices.
 %
@@ -79,7 +79,7 @@ end
 model_dynamic = str2func([M.fname,'_dynamic']);
 z = Y(find(lead_lag_incidence'));
 
-[d1,jacobian] = model_dynamic(z, exogenousvariables, params, steadystate,maximum_lag+1);
+[d1,jacobian] = model_dynamic(z, exogenousvariables, params, steadystate, oo.exo_steady_state, maximum_lag+1);
 
 res = zeros(periods*ny,1);
 
@@ -100,7 +100,7 @@ for iter = 1:options.simul.maxit
     i_cols = i_cols_A+(maximum_lag-1)*ny;
     m = 0;
     for it = (maximum_lag+1):(maximum_lag+periods)
-        [d1,jacobian] = model_dynamic(Y(i_cols), exogenousvariables, params, steadystate,it);
+        [d1,jacobian] = model_dynamic(Y(i_cols), exogenousvariables, params, steadystate, oo.exo_steady_state, it);
         if it == maximum_lag+periods && it == maximum_lag+1
             [r,c,v] = find(jacobian(:,i_cols_0));
             iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A0(c(:)),v(:)];
diff --git a/matlab/perfect-foresight-models/sim1_lbj.m b/matlab/perfect-foresight-models/sim1_lbj.m
index d4d8c24f2b..d5d3bd931d 100644
--- a/matlab/perfect-foresight-models/sim1_lbj.m
+++ b/matlab/perfect-foresight-models/sim1_lbj.m
@@ -1,4 +1,4 @@
-function [endogenousvariables, info] = sim1_lbj(endogenousvariables, exogenousvariables, steadystate, M, options)
+function [endogenousvariables, info] = sim1_lbj(endogenousvariables, exogenousvariables, steadystate, M, options, oo)
 
 % Performs deterministic simulations with lead or lag on one period using the historical LBJ algorithm
 %
@@ -15,7 +15,7 @@ function [endogenousvariables, info] = sim1_lbj(endogenousvariables, exogenousva
 % SPECIAL REQUIREMENTS
 %   None.
 
-% Copyright (C) 1996-2015 Dynare Team
+% Copyright (C) 1996-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -70,14 +70,14 @@ for iter = 1:options.simul.maxit
     end
     it_ = it_init;
     z = [endogenousvariables(iyp,it_-1) ; endogenousvariables(:,it_) ; endogenousvariables(iyf,it_+1)];
-    [d1, jacobian] = feval(dynamicmodel, z, exogenousvariables, M.params, steadystate, it_);
+    [d1, jacobian] = feval(dynamicmodel, z, exogenousvariables, M.params, steadystate, oo.exo_steady_state, it_);
     jacobian = [jacobian(:,iz), -d1];
     ic = [1:ny];
     icp = iyp;
     c (ic,:) = jacobian(:,is)\jacobian(:,isf1);
     for it_ = it_init+(1:options.periods-1)
         z = [endogenousvariables(iyp,it_-1) ; endogenousvariables(:,it_) ; endogenousvariables(iyf,it_+1)];
-        [d1, jacobian] = feval(dynamicmodel, z, exogenousvariables, M.params, steadystate, it_);
+        [d1, jacobian] = feval(dynamicmodel, z, exogenousvariables, M.params, steadystate, oo.exo_steady_state, it_);
         jacobian = [jacobian(:,iz), -d1];
         jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:);
         ic = ic + ny;
diff --git a/matlab/perfect-foresight-models/sim1_linear.m b/matlab/perfect-foresight-models/sim1_linear.m
index 6820c92bc2..dc057bfcf9 100644
--- a/matlab/perfect-foresight-models/sim1_linear.m
+++ b/matlab/perfect-foresight-models/sim1_linear.m
@@ -1,4 +1,4 @@
-function [endogenousvariables, info] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
+function [endogenousvariables, info] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options, oo)
 
 % Solves a linear approximation of a perfect foresight model using sparse matrix.
 %
@@ -38,7 +38,7 @@ function [endogenousvariables, info] = sim1_linear(endogenousvariables, exogenou
 % to center the variables around the deterministic steady state to solve the
 % perfect foresight model.
 
-% Copyright (C) 2015 Dynare Team
+% Copyright (C) 2015-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -111,7 +111,7 @@ dynamicmodel = str2func([M.fname,'_dynamic']);
 z = steadystate_y([ip; ic; in]);
 
 % Evaluate the Jacobian of the dynamic model at the deterministic steady state.
-[d1,jacobian] = dynamicmodel(z, transpose(steadystate_x), params, steadystate_y, 1);
+[d1,jacobian] = dynamicmodel(z, transpose(steadystate_x), params, steadystate_y, oo.exo_steady_state, 1);
 
 % Check that the dynamic model was evaluated at the steady state.
 if max(abs(d1))>1e-12
diff --git a/matlab/perfect-foresight-models/sim1_purely_backward.m b/matlab/perfect-foresight-models/sim1_purely_backward.m
index abc871a685..86121b54de 100644
--- a/matlab/perfect-foresight-models/sim1_purely_backward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_backward.m
@@ -1,8 +1,8 @@
-function [endogenousvariables, info] = sim1_purely_backward(endogenousvariables, exogenousvariables, steadystate, M, options)
+function [endogenousvariables, info] = sim1_purely_backward(endogenousvariables, exogenousvariables, steadystate, M, options, oo)
 
 % Performs deterministic simulation of a purely backward model
 
-% Copyright (C) 2012-2015 Dynare Team
+% Copyright (C) 2012-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -44,7 +44,7 @@ for it = 2:options.periods+1
     [tmp, check] = solve1(dynamicmodel, [yb1; yb], 1:M.endo_nbr, nyb+1:nyb+M.endo_nbr, ...
                           1, options.gstep, options.solve_tolf, options.solve_tolx, ...
                           options.simul.maxit, options.debug, exogenousvariables, ...
-                          M.params, steadystate, it+M.maximum_lag-1);
+                          M.params, steadystate, oo.exo_steady_state, it+M.maximum_lag-1);
     if check
         info.status = 0;
     end
diff --git a/matlab/perfect-foresight-models/sim1_purely_forward.m b/matlab/perfect-foresight-models/sim1_purely_forward.m
index ca2b2b5b27..81e7cb53ee 100644
--- a/matlab/perfect-foresight-models/sim1_purely_forward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_forward.m
@@ -1,7 +1,7 @@
-function [endogenousvariables, info] = sim1_purely_forward(endogenousvariables, exogenousvariables, steadystate, M, options)
+function [endogenousvariables, info] = sim1_purely_forward(endogenousvariables, exogenousvariables, steadystate, M, options, oo)
 % Performs deterministic simulation of a purely forward model
 
-% Copyright (C) 2012-2015 Dynare Team
+% Copyright (C) 2012-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,7 +36,7 @@ for it = options.periods:-1:1
                           1, options.gstep, options.solve_tolf, ...
                           options.solve_tolx, options.simul.maxit, ...
                           options.debug, exogenousvariables, M.params, steadystate, ...
-                          it+M.maximum_lag);
+                          oo.exo_steady_state, it+M.maximum_lag);
     if check
         info.status = 0;
     end
diff --git a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
index 52a2efe48c..7f635e9444 100644
--- a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
@@ -1,6 +1,6 @@
 function [endogenousvariables, info] = solve_stacked_linear_problem(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
 
-% Copyright (C) 2015 Dynare Team
+% Copyright (C) 2015-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -25,7 +25,7 @@ ic = find(M.lead_lag_incidence(2,:)');
 in = find(M.lead_lag_incidence(3,:)');
 
 % Evaluate the Jacobian of the dynamic model at the deterministic steady state.
-[d1,jacobian] = dynamicmodel(steadystate_y([ip; ic; in]), transpose(steadystate_x), M.params, steadystate_y, 1);
+[d1,jacobian] = dynamicmodel(steadystate_y([ip; ic; in]), transpose(steadystate_x), M.params, steadystate_y, steadystate_x, 1);
 
 % Check that the dynamic model was evaluated at the steady state.
 if max(abs(d1))>1e-12
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
index bff500b0d2..7545434ff9 100644
--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -1,4 +1,4 @@
-function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options)
+function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options, oo)
 % [endogenousvariables, info] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options);
 % Solves the perfect foresight model using dynare_solve
 %
@@ -48,14 +48,14 @@ if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity
     end
     [y, check] = dynare_solve(@perfect_foresight_mcp_problem,z(:),options, ...
                               dynamicmodel, y0, yT, ...
-                              exogenousvariables, M.params, steadystate, ...
+                              exogenousvariables, M.params, steadystate, oo.exo_steady_state, ...
                               M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
                               i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ...
                               M.NNZDerivatives(1),eq_index);
 else
     [y, check] = dynare_solve(@perfect_foresight_problem,z(:),options, ...
                               dynamicmodel, y0, yT, ...
-                              exogenousvariables, M.params, steadystate, ...
+                              exogenousvariables, M.params, steadystate, oo.exo_steady_state, ...
                               M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
                               i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ...
                               M.NNZDerivatives(1));
diff --git a/matlab/simul_backward_linear_model.m b/matlab/simul_backward_linear_model.m
index c4ed4fe264..5a3249ada0 100644
--- a/matlab/simul_backward_linear_model.m
+++ b/matlab/simul_backward_linear_model.m
@@ -34,7 +34,7 @@ function DynareOutput = simul_backward_linear_model(initial_conditions, sample_s
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 2012-2016 Dynare Team
+% Copyright (C) 2012-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -79,7 +79,7 @@ Y = DynareOutput.endo_simul;
 [cst,jacob] = model_dynamic(zeros(DynareModel.endo_nbr+ny1,1), ...
                             zeros(2,size(DynareOutput.exo_simul, 2)), ...
                             DynareModel.params, ...
-                            DynareOutput.steadystate,2);
+                            DynareOutput.steadystate, DynareOutput.exo_steady_state, 2);
 A0inv = inv(jacob(:,jdx));
 A1 = jacob(:,nonzeros(DynareModel.lead_lag_incidence(1,:)));
 B = jacob(:,end-number_of_shocks+1:end);
diff --git a/matlab/simul_backward_nonlinear_model.m b/matlab/simul_backward_nonlinear_model.m
index 99dd3932be..cb38a2af1c 100644
--- a/matlab/simul_backward_nonlinear_model.m
+++ b/matlab/simul_backward_nonlinear_model.m
@@ -34,7 +34,7 @@ function DynareOutput = simul_backward_nonlinear_model(initial_conditions, sampl
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 2012-2016 Dynare Team
+% Copyright (C) 2012-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -118,7 +118,7 @@ for it = 2:sample_size+1
                     DynareOptions.solve_tolf,DynareOptions.solve_tolx, ...
                     DynareOptions.simul.maxit,DynareOptions.debug, ...
                     DynareOutput.exo_simul, DynareModel.params, ...
-                    DynareOutput.steady_state, it);
+                    DynareOutput.steady_state, DynareOutput.exo_steady_state, it);
     Y(:,it) = z(jdx);
 end
 
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index 2772001431..3c3126d1ed 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -29,7 +29,7 @@ function [dr,info] = stochastic_solvers(dr,task,M_,options_,oo_)
 %   none.
 %  
 
-% Copyright (C) 1996-2016 Dynare Team
+% Copyright (C) 1996-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -110,7 +110,7 @@ if options_.order == 1
         jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd];
     else
         [junk,jacobia_] = feval([M_.fname '_dynamic'],z(iyr0),exo_simul, ...
-                            M_.params, dr.ys, it_);
+                            M_.params, dr.ys, oo_.exo_steady_state, it_);
     end;
 elseif options_.order == 2
     if (options_.bytecode)
@@ -120,7 +120,8 @@ elseif options_.order == 2
     else
         [junk,jacobia_,hessian1] = feval([M_.fname '_dynamic'],z(iyr0),...
                                          exo_simul, ...
-                                         M_.params, dr.ys, it_);
+                                         M_.params, dr.ys, oo_.exo_steady_state, ...
+                                         it_);
     end;
     if options_.use_dll
         % In USE_DLL mode, the hessian is in the 3-column sparse representation
diff --git a/matlab/thet2tau.m b/matlab/thet2tau.m
index 929193ca44..c34fd78e9d 100644
--- a/matlab/thet2tau.m
+++ b/matlab/thet2tau.m
@@ -1,7 +1,7 @@
 function tau = thet2tau(params, estim_params_, M_, oo_, indx, indexo, flagmoments,mf,nlags,useautocorr,iv)
 
 %
-% Copyright (C) 2011-2012 Dynare Team
+% Copyright (C) 2011-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -51,7 +51,7 @@ elseif flagmoments==-1
     [I,J]=find(M_.lead_lag_incidence');
     yy0=oo_.dr.ys(I);
     [residual, g1] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
-        M_.params, oo_.dr.ys, 1);
+        M_.params, oo_.dr.ys, oo_.exo_steady_state, 1);
     tau=[oo_.dr.ys(oo_.dr.order_var); g1(:)];
 
 else
diff --git a/tests/lmmcp/sw_newton.mod b/tests/lmmcp/sw_newton.mod
index 074fd82586..7fb5fc079e 100644
--- a/tests/lmmcp/sw_newton.mod
+++ b/tests/lmmcp/sw_newton.mod
@@ -69,11 +69,11 @@ end;
 
 simul(periods=1000);
 
-newton_solution_is_wrong = abs(evaluate_max_dynamic_residual(str2func('sw_newton_dynamic'), oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1000, size(oo_.endo_simul, 1), 1, M_.lead_lag_incidence))>options_.dynatol.f;
+newton_solution_is_wrong = abs(evaluate_max_dynamic_residual(str2func('sw_newton_dynamic'), oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, oo_.exo_steady_state, 1000, size(oo_.endo_simul, 1), 1, M_.lead_lag_incidence))>options_.dynatol.f;
 
 lmmcp = load('sw_lmmcp_results');
 
-lmmcp_solution_is_wrong = abs(evaluate_max_dynamic_residual(str2func('sw_newton_dynamic'), lmmcp.oo_.endo_simul, lmmcp.oo_.exo_simul, M_.params, oo_.steady_state, 1000, size(oo_.endo_simul, 1), 1, M_.lead_lag_incidence))>options_.dynatol.f;
+lmmcp_solution_is_wrong = abs(evaluate_max_dynamic_residual(str2func('sw_newton_dynamic'), lmmcp.oo_.endo_simul, lmmcp.oo_.exo_simul, M_.params, oo_.steady_state, oo_.exo_steady_state, 1000, size(oo_.endo_simul, 1), 1, M_.lead_lag_incidence))>options_.dynatol.f;
 
 solutions_are_different = max(max(abs(lmmcp.oo_.endo_simul-oo_.endo_simul)))>options_.dynatol.x;
 
-- 
GitLab