diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index 9e105274aaa95e066c4696e75f95bde7f334d44b..f730176db0823947f12d59183db1d7ad4433127b 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -192,9 +192,9 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
     warning(warning_old_state);
 end
 
-
+%If simulated paths are complex, take real part and recompute the residuals to check whether this is actually a solution
 if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_foresight_problem DLL
-    ny = size(oo_.endo_simul, 1)
+    ny = size(oo_.endo_simul, 1);
     if M_.maximum_lag > 0
         y0 = real(oo_.endo_simul(:, M_.maximum_lag));
     else
@@ -205,9 +205,20 @@ if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_fores
     else
         yT = NaN(ny, 1);
     end
-    yy = real(oo_.endo_simul(:,M_.maximum_lag+(1:periods)));
-
-    residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
+    if M_.maximum_lag~=0 && M_.maximum_lead~=0
+        yy = real(oo_.endo_simul(:,M_.maximum_lag+(1:periods)));
+        residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
+    else
+        %The perfect_foresight_problem MEX only works on models with lags and leads
+        i_cols = find(M_.lead_lag_incidence');
+        residuals=NaN(ny,periods);
+        yy=real(oo_.endo_simul);
+        for it = (M_.maximum_lag+1):(M_.maximum_lag+periods)
+            residuals(:,it) = feval([M_.fname '.dynamic'],yy(i_cols), oo_.exo_simul, M_.params, oo_.steady_state, it);
+            i_cols = i_cols + ny;
+        end
+        residuals=residuals(:);
+    end
 
     if max(abs(residuals))< options_.dynatol.f
         oo_.deterministic_simulation.status = true;