diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 712fd6ef1..ca98027a2 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -52,6 +52,7 @@ options_.lyapunov_complex_threshold = 1e-15;
 options_.solve_algo = 4;
 options_.solve_tolf = eps^(1/3);
 options_.solve_tolx = eps^(2/3);
+options_.solve_randomize_initial_guess = true;
 options_.trust_region_initial_step_bound_factor = 1;
 options_.dr_display_tol=1e-6;
 options_.minimal_workspace = false;
diff --git a/matlab/display_static_residuals.m b/matlab/display_static_residuals.m
index 0f14d533e..d6b62e81e 100644
--- a/matlab/display_static_residuals.m
+++ b/matlab/display_static_residuals.m
@@ -51,7 +51,7 @@ if any(imag(oo_.steady_state))
 end
 
 if options_.steadystate_flag
-    [oo_.steady_state,M_.params,info] = ...
+    [oo_.steady_state,M_.params] = ...
         evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,false);
 end
 
@@ -91,6 +91,8 @@ if nargout == 0
         first_eq = 1;
         last_eq = M_.orig_endo_nbr;
     end
+    disp_format_tags_real=sprintf('Equation number %%%uu: %%-%us: %%14.6f\n',length(num2str(M_.eq_nbr)),size(strvcat(tags(:,3)),2)+1);
+    disp_format_tags_complex=sprintf('Equation number %%%uu: %%-%us: real: %%14.6f, imaginary: %%g\n',length(num2str(M_.eq_nbr)),size(strvcat(tags(:,3)),2)+1);
     for i=first_eq:last_eq
         if abs(z(i)) < options_.solve_tolf/100
             tmp = 0;
@@ -104,11 +106,17 @@ if nargout == 0
         end
         if ~(non_zero && tmp == 0)
             if ~istag || length(ind) == 0
-                disp(['Equation number ' int2str(i) ' : ' num2str(tmp)])
+                if ~isreal(z)
+                    fprintf('Equation number %u: %g (imaginary part: %g)\n', i, real(tmp), imag(tmp));
+                else
+                    fprintf('Equation number %u: %g\n', i, tmp);
+                end
             else
-                t1 = tg( ind , 2 );
-                s = cell2mat(t1);
-                disp( ['Equation number ', int2str(i) ,' : ', num2str(tmp) ,' : ' s])
+                if ~isreal(z)
+                    fprintf(disp_format_tags_complex, i, tg{ind , 2}, real(tmp), imag(tmp) );
+                else
+                    fprintf(disp_format_tags_real, i, tg{ind , 2}, tmp);
+                end
             end
         end
     end
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 4fe2cc55e..42b1eb541 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -59,7 +59,6 @@ else
 end
 
 % checking initial values
-% TODO We should have an option to deactivate the randomization.
 if jacobian_flag
     [fvec, fjac] = feval(f, x, varargin{:});
     wrong_initial_guess_flag = false;
@@ -70,32 +69,38 @@ if jacobian_flag
             errorcode = -11;
             return;
         end
-        disp_verbose('Randomize initial guess...', options.verbosity)
-        % Let's try random numbers for the variables initialized with the default value.
-        wrong_initial_guess_flag = true;
-        % First try with positive numbers.
-        tentative_number = 0;
-        while wrong_initial_guess_flag && tentative_number<=in0*10
-            tentative_number = tentative_number+1;
-            x(idx) = rand(in0, 1)*10;
-            [fvec, fjac] = feval(f, x, varargin{:});
-            wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
-        end
-        % If all previous attempts failed, try with real numbers.
-        tentative_number = 0;
-        while wrong_initial_guess_flag && tentative_number<=in0*10
-            tentative_number = tentative_number+1;
-            x(idx) = randn(in0, 1)*10;
-            [fvec, fjac] = feval(f, x, varargin{:});
-            wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
-        end
-        % Last tentative, ff all previous attempts failed, try with negative numbers.
-        tentative_number = 0;
-        while wrong_initial_guess_flag && tentative_number<=in0*10
-            tentative_number = tentative_number+1;
-            x(idx) = -rand(in0, 1)*10;
-            [fvec, fjac] = feval(f, x, varargin{:});
-            wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
+        if options.solve_randomize_initial_guess
+            if any(~isreal(fvec)) || any(~isreal(fjac(:)))
+                disp_verbose('dynare_solve: starting value results in complex values. Randomize initial guess...', options.verbosity)
+            else
+                disp_verbose('dynare_solve: starting value results in nonfinite/NaN value. Randomize initial guess...', options.verbosity)
+            end
+            % Let's try random numbers for the variables initialized with the default value.
+            wrong_initial_guess_flag = true;
+            % First try with positive numbers.
+            tentative_number = 0;
+            while wrong_initial_guess_flag && tentative_number<=in0*10
+                tentative_number = tentative_number+1;
+                x(idx) = rand(in0, 1)*10;
+                [fvec, fjac] = feval(f, x, varargin{:});
+                wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)));
+            end
+            % If all previous attempts failed, try with real numbers.
+            tentative_number = 0;
+            while wrong_initial_guess_flag && tentative_number<=in0*10
+                tentative_number = tentative_number+1;
+                x(idx) = randn(in0, 1)*10;
+                [fvec, fjac] = feval(f, x, varargin{:});
+                wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)));
+            end
+            % Last tentative, ff all previous attempts failed, try with negative numbers.
+            tentative_number = 0;
+            while wrong_initial_guess_flag && tentative_number<=in0*10
+                tentative_number = tentative_number+1;
+                x(idx) = -rand(in0, 1)*10;
+                [fvec, fjac] = feval(f, x, varargin{:});
+                wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)));
+            end
         end
     end
 else
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 0f48675f3..49a33b2b9 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -293,22 +293,19 @@ elseif ~options.bytecode && ~options.block
                 options.mcppath.lb = lb;
                 options.mcppath.ub = ub;
             end
-            [ys,check,fvec] = dynare_solve(@static_mcp_problem,...
+            [ys,check,fvec, fjac, errorcode] = dynare_solve(@static_mcp_problem,...
                 ys_init,...
                 options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
                 options, exo_ss, params,...
                 M.endo_nbr, static_resid, static_g1, ...
                 M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, eq_index);
         else
-            [ys, check] = dynare_solve(@static_problem, ys_init, ...
+            [ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ...
                 options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
                 options, exo_ss, params, M.endo_nbr, static_resid, static_g1, ...
                 M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr);
         end
         if check && options.debug
-            [ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ...
-                                                              options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
-                                                              options, exo_ss, params, M.endo_nbr, static_resid, static_g1, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr);
             dprintf('Nonlinear solver routine returned errorcode=%i.', errorcode)
             skipline()
             [infrow,infcol]=find(isinf(fjac) | isnan(fjac));
@@ -480,10 +477,14 @@ if M.static_and_dynamic_models_differ
 end
 
 if ~isreal(ys)
-    info(1) = 21;
-    info(2) = sum(imag(ys).^2);
-    ys = real(ys);
+    if sum(imag(ys).^2) < 1e-7
+        ys=real(ys);
+    else
+        info(1) = 21;
+        info(2) = sum(imag(ys).^2);
+        ys = real(ys);
     return
+    end
 end
 
 if ~isempty(find(isnan(ys)))
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 4f57c869b..120d79908 100644
--- a/matlab/get_error_message.m
+++ b/matlab/get_error_message.m
@@ -67,9 +67,9 @@ switch info(1)
         message = 'The steadystate file did not compute the steady state';
     case 20
         if DynareOptions.linear
-            message = sprintf('Impossible to find the steady state (the sum of square residuals of the static equations is %5.4f). Either the model doesn''t have a steady state or there are an infinity of steady states Check whether your model is truly linear or whether there is a mistake in linearization.', info(2));
+            message = sprintf('Impossible to find the steady state (the sum of squared residuals of the static equations is %5.4f). Either the model doesn''t have a steady state or there are an infinity of steady states. Check whether your model is truly linear or whether there is a mistake in linearization.', info(2));
         else
-            message = sprintf('Impossible to find the steady state (the sum of square residuals of the static equations is %5.4f). Either the model doesn''t have a steady state, there are an infinity of steady states, or the guess values are too far from the solution', info(2));
+            message = sprintf('Impossible to find the steady state (the sum of squared residuals of the static equations is %5.4f). Either the model doesn''t have a steady state, there are an infinity of steady states, or the guess values are too far from the solution', info(2));
         end
     case 21
         message = sprintf('The steady state is complex (the sum of square residuals of imaginary parts of the steady state is %5.4f)', info(2));
diff --git a/matlab/steady.m b/matlab/steady.m
index bbbf070e8..255f4d7db 100644
--- a/matlab/steady.m
+++ b/matlab/steady.m
@@ -94,9 +94,18 @@ else
         end
     end
     if options_.debug
-        fprintf('\nThe steady state computation failed. It terminated with the following values:\n')
+        fprintf('\nsteady: The steady state computation failed. It terminated with the following values:\n')
+        if ~isreal(oo_.steady_state)
+            format_string=sprintf('%%-%us= %%g%%+gi\n',size(strvcat(M_.endo_names),2)+1);
+        else
+            format_string=sprintf('%%-%us= %%14.6f\n',size(strvcat(M_.endo_names),2)+1);
+        end
         for i=1:M_.orig_endo_nbr
-            fprintf('%s \t\t %g\n', M_.endo_names{i}, oo_.steady_state(i));
+            if ~isreal(oo_.steady_state)
+                fprintf(format_string, M_.endo_names{i}, real(oo_.steady_state(i)),imag(oo_.steady_state(i)));
+            else
+                fprintf(format_string, M_.endo_names{i}, oo_.steady_state(i));
+            end
         end
     end
     print_info(info,options_.noprint, options_);
diff --git a/preprocessor b/preprocessor
index d95c2f2bd..bdb5cdacf 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit d95c2f2bd81a4e8fc56552d19e62ac37b8d08a81
+Subproject commit bdb5cdacf22019c75db38dcdfe696f2978ab31f9