diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 72555ed859c44dc0d1741e6ef574244fd600a00a..d020df37027ca9209e8b54618ec00c2d2dd6df7a 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;
@@ -77,30 +105,6 @@ elseif options_.solve_algo == 1
     [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,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) ;
diff --git a/matlab/print_info.m b/matlab/print_info.m
index b406e964a004111787713a13f76ca84c7c827925..e197383815e03315005be447c276ad2ae0708c11 100644
--- a/matlab/print_info.m
+++ b/matlab/print_info.m
@@ -55,8 +55,8 @@ if ~noprint
         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 74111276f6468499888ea48c1f86f46792d75f02..83d3efb12e095acdb42a4e400cf0ae22c115a14d 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -58,9 +58,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 ;