diff --git a/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m b/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m
index 78d9fed418d182339acfa25b89e247e2346288a4..3a04ea8e1914dfc3c75cee2bad98c124381c8bc6 100644
--- a/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m
+++ b/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m
@@ -55,9 +55,9 @@ for it = maximum_lag+(1:T)
     z(jexog) = transpose(exo_simul(it,:));
     residuals(i_rows) = dynamicjacobian*z;
     if nargout == 2
-        if it == 2
+        if it == maximum_lag+1
             JJacobian(i_rows,i_cols_J1) = dynamicjacobian(:,i_cols_1);
-        elseif it == T + 1
+        elseif it == maximum_lag+T
             JJacobian(i_rows,i_cols_J(i_cols_T)) = dynamicjacobian(:,i_cols_T);
         else
             JJacobian(i_rows,i_cols_J) = dynamicjacobian(:,i_cols_j);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
index eb898322df33fce3da043103b0caabfa177c971b..000fff50f370fe9c7aa68e2f09f20e445cd07cd5 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m
@@ -74,7 +74,7 @@ i_rows = 1:ny;
 offset = 0;
 i_cols_J = i_cols;
 
-for it = 2:(T+1)
+for it = maximum_lag+(1:T)
     if nargout == 1
         res = dynamic_function(YY(i_cols),exo_simul, params, ...
                                steady_state,it);
@@ -83,15 +83,15 @@ for it = 2:(T+1)
         [res,jacobian] = dynamic_function(YY(i_cols),exo_simul, params, ...
                                           steady_state,it);
         residuals(i_rows) = res(eq_index);
-        if it == 2
+        if it == maximum_lag+1
             [rows,cols,vals] = find(jacobian(eq_index,i_cols_1));
             iJacobian{1} = [offset+rows, i_cols_J1(cols), vals];
-        elseif it == T + 1
+        elseif it == maximum_lag+T
             [rows,cols,vals] = find(jacobian(eq_index,i_cols_T));
             iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals];
         else
             [rows,cols,vals] = find(jacobian(eq_index,i_cols_j));
-            iJacobian{it-1} = [offset+rows, i_cols_J(cols), vals];
+            iJacobian{it-maximum_lag} = [offset+rows, i_cols_J(cols), vals];
             i_cols_J = i_cols_J + ny;
         end
         offset = offset + ny;
diff --git a/matlab/perfect-foresight-models/perfect_foresight_problem.m b/matlab/perfect-foresight-models/perfect_foresight_problem.m
index d27ec6eaf0115a39d8e160df60437d5a16363424..c6c6f77a71476242387e85b039612d2cf4098f4f 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_problem.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_problem.m
@@ -71,22 +71,22 @@ i_rows = 1:ny;
 i_cols_J = i_cols;
 offset = 0;
 
-for it = (maximum_lag+1):(T+1)
+for it = maximum_lag+(1:T)
     if nargout == 1
         residuals(i_rows) = dynamic_function(YY(i_cols),exo_simul, params, ...
                                              steady_state,it);
     elseif nargout == 2
         [residuals(i_rows),jacobian] = dynamic_function(YY(i_cols),exo_simul, params, ...
                                                         steady_state,it);
-        if it == 2
+        if it == maximum_lag+1
             [rows,cols,vals] = find(jacobian(:,i_cols_1));
             iJacobian{1} = [offset+rows, i_cols_J1(cols), vals];
-        elseif it == T + 1
+        elseif it == maximum_lag+T
             [rows,cols,vals] = find(jacobian(:,i_cols_T));
             iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals];
         else
             [rows,cols,vals] = find(jacobian(:,i_cols_j));
-            iJacobian{it-1} = [offset+rows, i_cols_J(cols), vals];
+            iJacobian{it-maximum_lag} = [offset+rows, i_cols_J(cols), vals];
             i_cols_J = i_cols_J + ny;
         end
         offset = offset + ny;
diff --git a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
index 014bb9acb44ac8000c679457c49c2a2c1e40c820..0f14b2367fc0328593da45aae1d3efa1253aad33 100644
--- a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
+++ b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
@@ -67,9 +67,9 @@ elseif (options.solve_algo == 11)
     end
 end
 
-y0 = endogenousvariables(:,1);
-yT = endogenousvariables(:,periods+2);
-z = endogenousvariables(:,2:periods+1);
+y0 = endogenousvariables(:,M.maximum_lag);
+yT = endogenousvariables(:,M.maximum_lag+periods+1);
+z = endogenousvariables(:,M.maximum_lag+(1:periods));
 illi = M.lead_lag_incidence';
 [i_cols, junk,i_cols_j] = find(illi(:));
 illi = illi(:,2:3);
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
index dd3bc6363d165f018e01e4bf5a7a8b78e4c056fb..7e6517be871bcd034245b0e6cebf7832d1fe3451 100644
--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -69,7 +69,7 @@ else
     check = 1;
 end
 
-endogenousvariables = [y0 reshape(y, M.endo_nbr, options.periods) yT];
+endogenousvariables(:, M.maximum_lag+(1:options.periods)) = reshape(y, M.endo_nbr, options.periods);
 
 if check
     info.status = false;
diff --git a/tests/deterministic_simulations/rbc_det.mod b/tests/deterministic_simulations/rbc_det.mod
index 71fdbd59f9e81d5b4ce2aba8f0a96c9d3c43e5fa..1454609c08eb13a29016231c5ea02251c5d52b31 100644
--- a/tests/deterministic_simulations/rbc_det.mod
+++ b/tests/deterministic_simulations/rbc_det.mod
@@ -17,7 +17,7 @@ sigma2  =   0;
 model;
 
   // Eq. n°1:
-  efficiency = rho*efficiency(-1) + EfficiencyInnovation;
+  efficiency = rho*efficiency(-1) + EfficiencyInnovation(-2); // Use a lag of two to test the maximum_lag logic
 
   // Eq. n°2:
   Efficiency = effstar*exp(efficiency);
diff --git a/tests/deterministic_simulations/rbc_det_stack_solve_algo_7.mod b/tests/deterministic_simulations/rbc_det_stack_solve_algo_7.mod
index a4c829ffb37d6102c2ff238544b60c863d9a43d1..baf60baf2901fc22635f98ace571e3583dcb666c 100644
--- a/tests/deterministic_simulations/rbc_det_stack_solve_algo_7.mod
+++ b/tests/deterministic_simulations/rbc_det_stack_solve_algo_7.mod
@@ -17,7 +17,7 @@ sigma2  =   0;
 model;
 
   // Eq. n°1:
-  efficiency = rho*efficiency(-1) + EfficiencyInnovation;
+  efficiency = rho*efficiency(-1) + EfficiencyInnovation(-2); // Use a lag of two to test the maximum_lag logic
 
   // Eq. n°2:
   Efficiency = effstar*exp(efficiency);