diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 8dfedc081..613b01112 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -42,6 +42,24 @@ nn = size(x,1);
 % checking initial values
 if jacobian_flag
     [fvec,fjac] = feval(func,x,varargin{:});
+    if any(any(isinf(fjac) | isnan(fjac)))
+        [infrow,infcol]=find(isinf(fjac) | isnan(fjac));
+        M=evalin('base','M_'); %get variable names from workspace
+        fprintf('\nSTEADY:  The Jacobian contains Inf or NaN. The problem arises from: \n\n')
+        for ii=1:length(infrow)
+            if infcol(ii)<=M.orig_endo_nbr
+                fprintf('STEADY:  Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',infrow(ii),M.endo_names(infcol(ii),:),M.endo_names(infcol(ii),:),x(infcol(ii)))
+            else %auxiliary vars
+                orig_var_index=M.aux_vars(1,infcol(ii)-M.orig_endo_nbr).orig_index;
+                fprintf('STEADY:  Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',infrow(ii),M.endo_names(orig_var_index,:),M.endo_names(orig_var_index,:),x(infcol(ii)))            
+            end
+        end
+        fprintf('\nSTEADY:  The problem most often occurs, because a variable with\n')
+        fprintf('STEADY:  exponent smaller than 1 has been initialized to 0. Taking the derivative\n')
+        fprintf('STEADY:  and evaluating it at the steady state then results in a division by 0.\n')
+        error('An element of the Jacobian is not finite or NaN') 
+    end
+
 else
     fvec = feval(func,x,varargin{:});
     fjac = zeros(nn,nn) ;