diff --git a/matlab/lnsrch1_wrapper_two_boundaries.m b/matlab/lnsrch1_wrapper_two_boundaries.m
index d1c5ecfcc..5b7bb5b39 100644
--- a/matlab/lnsrch1_wrapper_two_boundaries.m
+++ b/matlab/lnsrch1_wrapper_two_boundaries.m
@@ -47,5 +47,7 @@ function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, ...
 
 %reshape the input arguments of the dynamic function
 y(y_kmin+1:y_kmin+periods, y_index) = reshape(ya',length(y_index),periods)';
-[r, T, g1]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, y_size);
-ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*y_size, 1);
+ra = NaN(periods*y_size, 1);
+for it_ = y_kmin+(1:periods)
+    [ra((it_-y_kmin-1)*y_size+(1:y_size)), T, g1]=feval(fname, y, x, params, steady_state, T, it_, false);
+end
diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m
index f638975ad..aef8cf02c 100644
--- a/matlab/solve_two_boundaries.m
+++ b/matlab/solve_two_boundaries.m
@@ -80,11 +80,23 @@ ilu_setup.udiag = 0;
 max_resa=1e100;
 reduced = 0;
 while ~(cvg==1 || iter>maxit_)
-    [r, T, g1]=feval(fname, y, x, params, steady_state, T, periods, false, y_kmin, Blck_size);
+    r = NaN(Blck_size, periods);
+    g1a = spalloc(Blck_size*periods, Blck_size*periods, nze*periods);
+    for it_ = y_kmin+(1:periods)
+        [r(:, it_-y_kmin), T, g1]=feval(fname, y, x, params, steady_state, T, it_, false);
+        if periods == 1
+            g1a = g1(:, Blck_size+(1:Blck_size));
+        elseif it_ == y_kmin+1
+            g1a(1:Blck_size, 1:Blck_size*2) = g1(:, Blck_size+1:end);
+        elseif it_ == y_kmin+periods
+            g1a((periods-1)*Blck_size+1:end, (periods-2)*Blck_size+1:end) = g1(:, 1:2*Blck_size);
+        else
+            g1a((it_-y_kmin-1)*Blck_size+(1:Blck_size), (it_-y_kmin-2)*Blck_size+(1:3*Blck_size)) = g1;
+        end
+    end
     preconditioner = 2;
-    g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
-    ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)';
-    ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*Blck_size, 1);
+    ya = reshape(y(y_kmin+(1:periods),y_index)', 1, periods*Blck_size)';
+    ra = reshape(r, periods*Blck_size, 1);
     b=-ra+g1a*ya;
     [max_res, max_indx]=max(max(abs(r')));
     if ~isreal(r)
diff --git a/preprocessor b/preprocessor
index 85b900363..379be6cce 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 85b900363d58b42aa0d0dc4d5e3a2d02bf792448
+Subproject commit 379be6ccef55dc38af04edb881bf6aa57bdecd65