diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index da0e58921ada2e5c4372ffdfd133ca70cfff2669..5e9b7e89a487e070716e7f113b4b1106a9c1e066 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -35,7 +35,35 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
 
 global options_
 
+tolf = options_.solve_tolf ;
 info = 0;
+
+% checking initial values
+if jacobian_flag
+    [fvec,fjac] = feval(func,x,varargin{:});
+else
+    fvec = feval(func,x,varargin{:});
+    fjac = zeros(nn,nn) ;
+end
+
+i = find(~isfinite(fvec));
+
+if ~isempty(i)
+    disp(['STEADY:  numerical initial values or parameters incompatible with the following' ...
+          ' equations'])
+    disp(i')
+    disp('Please check for example')
+    disp('   i) if all parameters occurring in these equations are defined')
+    disp('  ii) that no division by an endogenous variable initialized to 0 occurs') 
+    info = 1;
+    x = NaN;
+    return;
+end
+
+if max(abs(fvec)) < tolf
+    return ;
+end
+
 if options_.solve_algo == 0
     if ~exist('OCTAVE_VERSION')
         if ~user_has_matlab_license('optimization_toolbox')
@@ -60,7 +88,7 @@ if options_.solve_algo == 0
         func = @(x) func2(x, varargin{:});
         % The Octave version of fsolve does not converge when it starts from the solution
         fvec = feval(func,x);
-        if max(abs(fvec)) >= options_.solve_tolf
+        if max(abs(fvec)) >= tolf
             [x,fval,exitval,output] = fsolve(func,x,options);
         else
             exitval = 3;
@@ -75,34 +103,10 @@ if options_.solve_algo == 0
 elseif options_.solve_algo == 1
     nn = size(x,1);
     [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,options_.gstep, ...
-                    options_.solve_tolf,options_.solve_tolx, ...
+                    tolf,options_.solve_tolx, ...
                     options_.solve_maxit,options_.debug,varargin{:});
 elseif options_.solve_algo == 2 || options_.solve_algo == 4
     nn = size(x,1) ;
-    tolf = options_.solve_tolf ;
-
-    if jacobian_flag
-        [fvec,fjac] = feval(func,x,varargin{:});
-    else
-        fvec = feval(func,x,varargin{:});
-        fjac = zeros(nn,nn) ;
-    end
-
-    i = find(~isfinite(fvec));
-    
-    if ~isempty(i)
-        disp(['STEADY:  numerical initial values incompatible with the following' ...
-              ' equations'])
-        disp(i')
-        disp('Please check for example')
-        disp('   i) if all parameters occurring in these equations are defined')
-        disp('  ii) that no division by an endogenous variable initialized to 0 occurs') 
-        error('exiting ...')
-    end
-    
-    if max(abs(fvec)) < tolf
-        return ;
-    end
 
     if ~jacobian_flag
         fjac = zeros(nn,nn) ;
@@ -129,7 +133,7 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
         end
         [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, ...
                         bad_cond_flag, options_.gstep, ...
-                        options_.solve_tolf,options_.solve_tolx, ...
+                        tolf,options_.solve_tolx, ...
                         options_.solve_maxit,options_.debug,varargin{:});
         if info
             return
@@ -138,7 +142,7 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
     fvec = feval(func,x,varargin{:});
     if max(abs(fvec)) > tolf
         [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, bad_cond_flag, ...
-                        options_.gstep, options_.solve_tolf,options_.solve_tolx, ...
+                        options_.gstep, tolf,options_.solve_tolx, ...
                         options_.solve_maxit,options_.debug,varargin{:});
     end
 elseif options_.solve_algo == 3
diff --git a/matlab/print_info.m b/matlab/print_info.m
index 15a793614fe22d181594d8ed8cc4380b175e69c6..a1245ce5e0479fbacbac4ec4c6134257937b12ad 100644
--- a/matlab/print_info.m
+++ b/matlab/print_info.m
@@ -64,12 +64,12 @@ if ~noprint
           error(['The Jacobian contains NaNs'])
         end
 
-        case 19
+      case 19
         error('The steadystate file did not compute the steady state')
       case 20
         error(['Impossible to find the steady state. Either the model' ...
-               ' doesn''t have a unique steady state of the guess values' ...
-               ' are too far from the solution'])
+               ' doesn''t have a steady state or there is an infinity of steady states' ...
+               ' or the guess values are too far from the solution'])
       case 21
         error('The steady state is complex')
       case 22
diff --git a/matlab/solve1.m b/matlab/solve1.m
index 97cb5e24e3fab5514f7c29f3010d5aa0c5f06bd7..52eb8460e6193b1b84179b54b875290216068113 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -59,9 +59,11 @@ fvec = fvec(j1);
 i = find(~isfinite(fvec));
 
 if ~isempty(i)
-    disp(['STEADY:  numerical initial values incompatible with the following' ...
-          ' equations'])
+    disp(['SOLVE1: during the resolution of the non-linear system, the evaluation of the following ' ...
+          'equation(s) resulted in a non-finite number:'])
     disp(j1(i)')
+    check = 1;
+    return;
 end
 
 f = 0.5*fvec'*fvec ;