diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index ebb6c8700c35132848be9333186e1e2fb7505338..3fc811add33186d9c458a2cc8ad49733dedd72c6 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -207,12 +207,12 @@ while (t <= sample_size)
             exo_simul(2,:) = exo_simul_(M_.maximum_lag+t,:) + ...
                 shocks((t-2)*replic_nbr+k,:);
             initial_conditions = results{k}(:,t-1);
-            results{k}(:,t) = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
-                                                 exo_simul,ep.init,initial_conditions,...
-                                                 maximum_lag,maximum_lead,steady_state, ...
-                                                 ep.verbosity,bytecode_flag,ep.stochastic.order,...
-                                                 M_.params,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
-                                                 options_.lmmcp,options_,oo_);
+            [results{k}(:,t), info_convergence] = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
+                                                              exo_simul,ep.init,initial_conditions,...
+                                                              maximum_lag,maximum_lead,steady_state, ...
+                                                              ep.verbosity,bytecode_flag,ep.stochastic.order,...
+                                                              M_.params,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
+                                                              options_.lmmcp,options_,oo_);
         end
     else
         for k = 1:replic_nbr
@@ -223,16 +223,21 @@ while (t <= sample_size)
             exo_simul(maximum_lag+1,:) = ...
                 exo_simul_(maximum_lag+t,:) + shocks((t-2)*replic_nbr+k,:);
             initial_conditions = results{k}(:,t-1);
-            results{k}(:,t) = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
-                                                 exo_simul,ep.init,initial_conditions,...
-                                                 maximum_lag,maximum_lead,steady_state, ...
-                                                 ep.verbosity,bytecode_flag,ep.stochastic.order,...
-                                                 M_,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
-                                                 options_.lmmcp,options_,oo_);
+            [results{k}(:,t), info_convergence] = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
+                                                              exo_simul,ep.init,initial_conditions,...
+                                                              maximum_lag,maximum_lead,steady_state, ...
+                                                              ep.verbosity,bytecode_flag,ep.stochastic.order,...
+                                                              M_,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
+                                                              options_.lmmcp,options_,oo_);
+        end
+    end
+    if verbosity
+        if info_convergence
+            disp(['Time: ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
+        else
+            disp(['Time: ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
         end
     end
-            
-    
 end% (while) loop over t
 
 dyn_waitbar_close(hh);
@@ -263,13 +268,12 @@ end
  assignin('base', 'Simulated_time_series', ts);
  
  
-function y = extended_path_core(periods,endo_nbr,exo_nbr,positive_var_indx, ...
+function [y, info_convergence] = extended_path_core(periods,endo_nbr,exo_nbr,positive_var_indx, ...
                                 exo_simul,init,initial_conditions,...
                                 maximum_lag,maximum_lead,steady_state, ...
                                 verbosity,bytecode_flag,order,M,pfm,algo,solve_algo,stack_solve_algo,...
                                 olmmcp,options,oo)
     
-    
 ep = options.ep;
 if init% Compute first order solution (Perturbation)...
     endo_simul = simult_(initial_conditions,oo.dr,exo_simul(2:end,:),1);
@@ -321,23 +325,16 @@ if flag
     else
         switch(algo)
           case 0
-            [flag,tmp] = ...
+            [flag,endo_simul] = ...
                 solve_stochastic_perfect_foresight_model(endo_simul,exo_simul,pfm,ep.stochastic.quadrature.nodes,ep.stochastic.order);
           case 1
-            [flag,tmp] = ...
+            [flag,endo_simul] = ...
                 solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,options_,pfm,ep.stochastic.order);
         end
-        endo_simul = tmp;
+        tmp.endo_simul = endo_simul;
         info_convergence = ~flag;
     end
 end
-if verbosity
-    if info_convergence
-        disp(['Time: ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
-    else
-        disp(['Time: ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
-    end
-end
 if info_convergence
     y = tmp.endo_simul(:,2);
 else