diff --git a/matlab/extended_path.m b/matlab/extended_path.m
index 462997e9a..7887a1974 100644
--- a/matlab/extended_path.m
+++ b/matlab/extended_path.m
@@ -60,6 +60,7 @@ if init
     oldopt = options_;
     options_.order = 1;
     [dr,info]=resol(oo_.steady_state,0);
+    oo_.dr = dr;
     options_ = oldopt;
     if init==2
         lambda = .8;
@@ -96,14 +97,18 @@ for t=1:sample_size
             oo_.endo_simul = initial_path(:,1:end-1)*lambda + oo_.endo_simul*(1-lambda);  
         end
     end
-    info = perfect_foresight_simulation;
+    if init  
+        info = perfect_foresight_simulation(dr,oo_.steady_state);
+    else
+        info = perfect_foresight_simulation;
+    end
     time = info.time;
     if verbose
         t
         info
     end
     if ~info.convergence
-        info = homotopic_steps(tdx,positive_var_indx,shocks,norme,.5);
+        info = homotopic_steps(tdx,positive_var_indx,shocks,norme,.5,init);
         if verbose
             norme
             info
@@ -122,7 +127,7 @@ end
  
  
 
-function info = homotopic_steps(tdx,positive_var_indx,shocks,init_weight,step)
+function info = homotopic_steps(tdx,positive_var_indx,shocks,init_weight,step,init)
     global oo_
     weight   = init_weight;
     verbose  = 1;
@@ -134,7 +139,11 @@ function info = homotopic_steps(tdx,positive_var_indx,shocks,init_weight,step)
         old_weight = weight;
         weight = weight+step;
         oo_.exo_simul(tdx,positive_var_indx) = weight*shocks+(1-weight);
-        info = perfect_foresight_simulation;
+        if init  
+            info = perfect_foresight_simulation(oo_.dr,oo_.steady_state);
+        else
+            info = perfect_foresight_simulation;
+        end
         time = time+info.time;
         if verbose
             [iter,step]
@@ -157,11 +166,15 @@ function info = homotopic_steps(tdx,positive_var_indx,shocks,init_weight,step)
     end
     if reduce_step
         step=step/1.5;
-        info = homotopic_steps(tdx,positive_var_indx,shocks,old_weight,step);
+        info = homotopic_steps(tdx,positive_var_indx,shocks,old_weight,step,init);
         time = time+info.time;
     elseif weight<1 && iter<100
         oo_.exo_simul(tdx,positive_var_indx) = shocks;
-        info = perfect_foresight_simulation;
+        if init  
+            info = perfect_foresight_simulation(oo_.dr,oo_.steady_state);
+        else
+            info = perfect_foresight_simulation;
+        end
         info.time = info.time+time;
     else
         info.time = time;
diff --git a/matlab/make_y_.m b/matlab/make_y_.m
index 1ce426e8e..37c731685 100644
--- a/matlab/make_y_.m
+++ b/matlab/make_y_.m
@@ -48,7 +48,7 @@ elseif size(oo_.endo_simul,2) < M_.maximum_lag+M_.maximum_lead+options_.periods
       case 0
         oo_.endo_simul = [oo_.endo_simul ...
                           oo_.steady_state*ones(1,M_.maximum_lag+options_.periods+M_.maximum_lead-size(oo_.endo_simul,2),1)];
-      case 1% A linear approximation is used to initiate the solution.
+      case 1% A linear approximation is used to initialize the solution.
         oldopt = options_;
         options_.order = 1;
         dr = oo_.dr;
diff --git a/matlab/perfect_foresight_simulation.m b/matlab/perfect_foresight_simulation.m
index 214e742f0..ea897ada6 100644
--- a/matlab/perfect_foresight_simulation.m
+++ b/matlab/perfect_foresight_simulation.m
@@ -70,8 +70,14 @@ else
     end
 end
 
-if compute_linear_solution
+if ~isstruct(compute_linear_solution) && compute_linear_solution 
     [dr,info]=resol(steady_state,0);
+elseif isstruct(compute_linear_solution)
+    dr = compute_linear_solution;
+    compute_linear_solution = 1;
+end
+
+if compute_linear_solution
     ghx = dr.ghx(end-dr.nfwrd+1:end,:);
 end
 
@@ -104,7 +110,7 @@ for iter = 1:options_.maxit_
     ic = 1:ny;
     icp = iyp;
     c(ic,:) = jacobian(:,is)\jacobian(:,isf1) ; 
-    for it_ = it_init+(1:periods-1)
+    for it_ = it_init+(1:periods-1-(options_.terminal_condition==2))
         z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)]; 
         [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_); 
         jacobian = [jacobian(:,iz) , -d1];
@@ -120,6 +126,7 @@ for iter = 1:options_.maxit_
             ic = ic + ny;
             c(ic,nrc) = s\c(ic,nrc);
         else% Terminal condition is Y_{T+1}-Y^{\star} = TransitionMatrix*(Y_{T}-Y^{\star})
+            it_ = it_+1;
             z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1) ] ;
             [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_);
             jacobian = [jacobian(:,iz) -d1];
@@ -127,16 +134,12 @@ for iter = 1:options_.maxit_
             ic = ic + ny;
             icp = icp + ny;
             s = jacobian(:,is);
-            iyp
-            nyp
-            isf
-            iyp-nyp
             s(:,iyp-nyp) = s(:,iyp-nyp)+jacobian(:,isf)*ghx;
             c (ic,:) = s\jacobian(:,isf1);
         end
         c = bksup0(c,ny,nrc,iyf,icf,periods);
         c = reshape(c,ny,periods+1);
-        endo_simul(:,it_init+(0:periods)) = endo_simul(:,it_init+(0:periods))+options_.slowc*c;
+        oo_.endo_simul(:,it_init+(0:periods)) = oo_.endo_simul(:,it_init+(0:periods))+options_.slowc*c;
     else% Terminal condition is Y_{T}=Y^{\star}
         c = bksup0(c,ny,nrc,iyf,icf,periods);
         c = reshape(c,ny,periods);