diff --git a/matlab/perfect-foresight-models/evaluate_max_dynamic_residual.m b/matlab/perfect-foresight-models/evaluate_max_dynamic_residual.m
deleted file mode 100644
index 22befabb6d132f8af678d832e5dbd8ba76a8cd39..0000000000000000000000000000000000000000
--- a/matlab/perfect-foresight-models/evaluate_max_dynamic_residual.m
+++ /dev/null
@@ -1,33 +0,0 @@
-function err = evaluate_max_dynamic_residual(model_dynamic, Y, exogenous_variables, params, steady_state, periods, ny, max_lag, lead_lag_incidence)
-
-% Copyright © 2013-2017 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/>.
-
-i_rows = 1:ny;
-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);
-    i_rows = i_rows + ny;
-    i_cols = i_cols + ny;
-    r = max(abs(d));
-    if r>err
-        err = r;
-    end
-end
\ No newline at end of file
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 6fae2f151e175efe99bbc3d0b39a0fdee4b6a06e..f96cb8832b921bec0b2d1e90d43060924a534e85 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -142,7 +142,7 @@ endogenousvariables(:, M_.maximum_lag+(1:periods)) = reshape(y, ny, periods);
 
 if options_.endogenous_terminal_period
     periods = options_.periods;
-    err = evaluate_max_dynamic_residual(str2func([M_.fname,'.dynamic']), endogenousvariables, exogenousvariables, M_.params, steadystate, periods, ny, M_.maximum_endo_lag, M_.lead_lag_incidence);
+    err = evaluate_max_dynamic_residual(str2func([M_.fname,'.sparse.dynamic_resid']), endogenousvariables, exogenousvariables, M_.params, steadystate, periods, M_.maximum_lag);
 end
 
 if stop
@@ -357,3 +357,14 @@ if rank_jacob < size(jacob,1)
         end
     end
 end
+
+
+function err = evaluate_max_dynamic_residual(dynamic_resid, endogenousvariables, exogenousvariables, params, steady_state, periods, max_lag)
+
+err = 0;
+
+for it = max_lag+(1:periods)
+    d = dynamic_resid(endogenousvariables(:, it+(-1:1)), exogenousvariables(it, :), params, steady_state);
+    r = max(abs(d));
+    err = max(err, r);
+end
diff --git a/tests/lmmcp/sw_newton.mod b/tests/lmmcp/sw_newton.mod
index 6162c796771c8744e136edd7141291eba068e05e..c450791c34a55ec9a5c77a6eef6235458c7d02b0 100644
--- a/tests/lmmcp/sw_newton.mod
+++ b/tests/lmmcp/sw_newton.mod
@@ -70,11 +70,15 @@ end;
 perfect_foresight_setup(periods=1000);
 perfect_foresight_solver;
 
-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;
 
-lmmcp = load(['sw_lmmcp' filesep 'Output' filesep 'sw_lmmcp_results']);
+verbatim;
+
+newton_resid = perfect_foresight_problem(reshape(oo_.endo_simul(:, 2:end-1),[],1), oo_.endo_simul(:, 1), oo_.endo_simul(:, end), oo_.exo_simul, M_.params, oo_.steady_state, options_.periods, M_, options_);
+newton_solution_is_wrong = max(abs(newton_resid)) > 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, 1000, size(oo_.endo_simul, 1), 1, M_.lead_lag_incidence))>options_.dynatol.f;
+lmmcp = load(['sw_lmmcp' filesep 'Output' filesep 'sw_lmmcp_results']);
+lmmcp_resid = perfect_foresight_problem(reshape(lmmcp.oo_.endo_simul(:, 2:end-1),[],1), lmmcp.oo_.endo_simul(:, 1), lmmcp.oo_.endo_simul(:, end), lmmcp.oo_.exo_simul, M_.params, lmmcp.oo_.steady_state, options_.periods, M_, options_);
+lmmcp_solution_is_wrong = max(abs(lmmcp_resid)) > options_.dynatol.f;
 
 solutions_are_different = max(max(abs(lmmcp.oo_.endo_simul-oo_.endo_simul)))>options_.dynatol.x;
 
@@ -89,3 +93,5 @@ end
 if lmmcp_solution_is_wrong
    error('Failed to solve SW with ZLB (using LMMCP algorithm on stacked model)!')
 end
+
+end; // verbatim block