diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index c23f562adc7aab1c8f4fbc369dca0281e86df8b2..0f48675f35d506ab0f192fff53e162dc99ccbb81 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -56,6 +56,15 @@ if isoctave && options.solve_algo == 11
     error(['STEADY: you can''t use solve_algo = %u under Octave'],options.solve_algo)
 end
 
+% To ensure that the z and zx matrices constructed by repmat and passed to bytecode
+% are of the right size.
+if size(ys_init, 2) > 1
+    error('ys_init must be a column-vector')
+end
+if size(exo_ss, 2) > 1
+    error('exo_ss must be a column-vector')
+end
+
 info = 0;
 check = 0;
 
@@ -458,7 +467,7 @@ if M.static_and_dynamic_models_differ
     % computed on the *static* one
     if options.bytecode
         z = repmat(ys,1,M.maximum_lead + M.maximum_lag + 1);
-        zx = repmat([exo_ss'], M.maximum_lead + M.maximum_lag + 1, 1);
+        zx = repmat(exo_ss', M.maximum_lead + M.maximum_lag + 1, 1);
         r = bytecode('dynamic','evaluate', z, zx, params, ys, 1);
     else
         r = feval([M.fname '.sparse.dynamic_resid'], repmat(ys, 3, 1), exo_ss, params, ys);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index d412afb30122659873708df3efbefbb9166657db..97f53848fd5e52a42a340f7edff1e4b0175c6c87 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -108,7 +108,7 @@ end
 if ~options_.simul.endval_steady && ~isempty(ys0_)
     terminal_condition_is_a_steady_state = true;
     for j = lastperiods
-        endval_resid = evaluate_static_model(oo_.endo_simul(:,j), oo_.exo_simul(j,:), M_.params, M_, options_);
+        endval_resid = evaluate_static_model(oo_.endo_simul(:,j), oo_.exo_simul(j,:)', M_.params, M_, options_);
         if norm(endval_resid, 'Inf') > options_.simul.steady_tolf
             terminal_condition_is_a_steady_state = false;
             break
@@ -160,7 +160,7 @@ function local_success = create_scenario(share)
         % Effectively compute the terminal steady state
         for j = lastperiods
             % First use the terminal steady of the previous homotopy iteration as guess value (or the contents of the endval block if this is the first iteration)
-            [oo_.endo_simul(:, j), ~, info] = evaluate_steady_state(oo_.endo_simul(:, j), oo_.exo_simul(j, :), M_, options_, true);
+            [oo_.endo_simul(:, j), ~, info] = evaluate_steady_state(oo_.endo_simul(:, j), oo_.exo_simul(j, :)', M_, options_, true);
             if info(1)
                 % If this fails, then try again using the initial steady state as guess value
                 if isempty(ys0_)
@@ -168,7 +168,7 @@ function local_success = create_scenario(share)
                 else
                     guess_value = ys0_;
                 end
-                [oo_.endo_simul(:, j), ~, info] = evaluate_steady_state(guess_value, oo_.exo_simul(j, :), M_, options_, true);
+                [oo_.endo_simul(:, j), ~, info] = evaluate_steady_state(guess_value, oo_.exo_simul(j, :)', M_, options_, true);
                 if info(1)
                     % If this fails again, give up and restore last periods in oo_.endo_simul
                     oo_.endo_simul(:, lastperiods) = saved_ss;