diff --git a/matlab/display_problematic_vars_Jacobian.m b/matlab/display_problematic_vars_Jacobian.m
new file mode 100644
index 0000000000000000000000000000000000000000..0cb42b581a4edb613dbe911a1a94b0fe800f352e
--- /dev/null
+++ b/matlab/display_problematic_vars_Jacobian.m
@@ -0,0 +1,73 @@
+function []=display_problematic_vars_Jacobian(problemrow,problemcol,M_,x,type,caller_string)
+% []=display_problematic_vars_Jacobian(problemrow,problemcol,M_,ys,caller_string)
+% print the equation numbers and variables associated with problematic entries 
+% of the Jacobian 
+%
+% INPUTS
+%   problemrow      [vector] rows associated with problematic entries
+%   problemcol      [vector] columns associated with problematic entries
+%   M_              [matlab structure] Definition of the model.           
+%   x               [vector] point at which the Jacobian was evaluated
+%   type            [string] 'static' or 'dynamic' depending on the type of
+%                               Jacobian
+%   caller_string   [string] contains name of calling function for printing 
+%    
+% OUTPUTS
+%   none.
+%  
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+if nargin<6
+    caller_string='';
+end
+if strcmp(type,'dynamic')
+    for ii=1:length(problemrow)
+        [var_row,var_index]=find(M_.lead_lag_incidence==problemcol(ii));
+        if var_row==2
+            type_string='';
+        elseif var_row==1
+            type_string='lag of';
+        elseif var_row==3;
+            type_string='lead of';
+        end
+        if var_index<=M_.orig_endo_nbr
+            fprintf('Derivative of Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',problemrow(ii),type_string,deblank(M_.endo_names(var_index,:)),deblank(M_.endo_names(var_index,:)),x(var_index))
+        else %auxiliary vars
+            orig_var_index=M_.aux_vars(1,var_index-M_.orig_endo_nbr).orig_index;
+            fprintf('Derivative of Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',problemrow(ii),type_string,deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),x(orig_var_index))            
+        end    
+    end
+    fprintf('\n%s  The problem most often occurs, because a variable with\n',caller_string)
+    fprintf('%s  exponent smaller than 1 has been initialized to 0. Taking the derivative\n',caller_string)
+    fprintf('%s  and evaluating it at the steady state then results in a division by 0.\n',caller_string)
+elseif strcmp(type,'static')
+    for ii=1:length(problemrow)
+        if problemcol(ii)<=M_.orig_endo_nbr
+            fprintf('Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',problemrow(ii),deblank(M_.endo_names(problemcol(ii),:)),deblank(M_.endo_names(problemcol(ii),:)),x(problemcol(ii)))
+        else %auxiliary vars
+            orig_var_index=M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).orig_index;
+            fprintf('Derivative of Equation %d with respect to Variable %s  (initial value of %s: %g) \n',problemrow(ii),deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),x(problemcol(ii)))            
+        end
+    end
+    fprintf('\n%s  The problem most often occurs, because a variable with\n',caller_string)
+    fprintf('%s exponent smaller than 1 has been initialized to 0. Taking the derivative\n',caller_string)
+    fprintf('%s and evaluating it at the steady state then results in a division by 0.\n',caller_string)
+else
+    error('Unknown Type')    
+end
\ No newline at end of file
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index b1ab7aefd9cb1232a448131c9caf7e53b99b5e13..744d174e2cf019ed30bfc5b32a2e3d2fd022936e 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -46,20 +46,9 @@ if jacobian_flag
         [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),deblank(M.endo_names(infcol(ii),:)),deblank(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),deblank(M.endo_names(orig_var_index,:)),deblank(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')
+        display_problematic_vars_Jacobian(infrow,infcol,M,x,'static','STEADY: ')
         error('An element of the Jacobian is not finite or NaN') 
     end
-
 else
     fvec = feval(func,x,varargin{:});
     fjac = zeros(nn,nn) ;
diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m
index e49eb29ab979ac09c726e6e2f60f42a509b38d55..197953fe1eb267876d250397a42796df76dd0e5b 100644
--- a/matlab/model_diagnostics.m
+++ b/matlab/model_diagnostics.m
@@ -47,7 +47,7 @@ problem_dummy=0;
 k = find(lead_lag_incidence(maximum_endo_lag+1,:)==0);
 if ~isempty(k)
     problem_dummy=1;
-    disp(['The following endogenous variables aren''t present at ' ...
+    disp(['MODEL_DIAGNOSTICS: The following endogenous variables aren''t present at ' ...
           'the current period in the model:'])
     for i=1:length(k)
         disp(endo_names(k(i),:))
@@ -64,24 +64,25 @@ if M.exo_nbr == 0
 end
 
 % check if ys is steady state
+options.debug=1; %locally set debug option to 1
 [dr.ys,params,check1]=evaluate_steady_state(oo.steady_state,M,options,oo,1);
 
 % testing for problem
 if check1(1)
     problem_dummy=1;
-    disp('model diagnostic can''t obtain the steady state')
+    disp('MODEL_DIAGNOSTICS: The steady state cannot be computed')
     if any(isnan(dr.ys))
-        disp(['model diagnostic obtains a steady state with NaNs'])
+        disp(['MODEL_DIAGNOSTICS: Steady state contains NaNs'])
     end
     if any(isinf(dr.ys))
-        disp(['model diagnostic obtains a steady state with Inf'])
+        disp(['MODEL_DIAGNOSTICS: Steady state contains Inf'])
     end
     return;
 end
 
 if ~isreal(dr.ys)
     problem_dummy=1;
-    disp(['model diagnostic obtains a steady state with complex ' ...
+    disp(['MODEL_DIAGNOSTICS: Steady state contains complex ' ...
           'numbers'])
     return
 end
@@ -110,13 +111,23 @@ for b=1:nb
     else
         [res,jacob]=feval([M.fname '_static'],dr.ys,exo,M.params);
     end
-    rank_jacob = rank(jacob);
+    if any(any(isinf(jacob) | isnan(jacob)))
+        problem_dummy=1;
+        [infrow,infcol]=find(isinf(jacob) | isnan(jacob));
+        fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the static model contains Inf or NaN. The problem arises from: \n\n')
+        display_problematic_vars_Jacobian(infrow,infcol,M,dr.ys,'static','MODEL_DIAGNOSTICS: ')
+    end
+    try
+        rank_jacob = rank(jacob); %can sometimes fail
+    catch
+        rank_jacob=size(jacob,1);
+    end
     if rank_jacob < size(jacob,1)
         problem_dummy=1;
         singularity_problem = 1;
-        disp(['model_diagnostic: the Jacobian of the static model is ' ...
+        disp(['MODEL_DIAGNOSTICS:  The Jacobian of the static model is ' ...
               'singular'])
-        disp(['there is ' num2str(endo_nbr-rank_jacob) ...
+        disp(['MODEL_DIAGNOSTICS:  there is ' num2str(endo_nbr-rank_jacob) ...
               ' colinear relationships between the variables and the equations'])
         ncol = null(jacob);
         n_rel = size(ncol,2);
@@ -150,13 +161,60 @@ for b=1:nb
         end
     end
 end
+    
 if singularity_problem 
-    fprintf('The presence of a singularity problem typically indicates that there is one\n')
-    fprintf('redundant equation entered in the model block, while another non-redundant equation\n')
-    fprintf('is missing. The problem often derives from Walras Law.\n')
+    fprintf('MODEL_DIAGNOSTICS:  The presence of a singularity problem typically indicates that there is one\n')
+    fprintf('MODEL_DIAGNOSTICS:  redundant equation entered in the model block, while another non-redundant equation\n')
+    fprintf('MODEL_DIAGNOSTICS:  is missing. The problem often derives from Walras Law.\n')
+end
+
+%%check dynamic Jacobian
+klen = M.maximum_lag + M.maximum_lead + 1;
+exo_simul = [repmat(oo.exo_steady_state',klen,1) repmat(oo.exo_det_steady_state',klen,1)];
+iyv = M.lead_lag_incidence';
+iyv = iyv(:);
+iyr0 = find(iyv) ;
+it_ = M.maximum_lag + 1;
+z = repmat(dr.ys,1,klen);
+
+if options.order == 1
+    if (options.bytecode)
+        [chck, junk, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
+                                        M.params, dr.ys, 1);
+        jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd];
+    else
+        [junk,jacobia_] = feval([M.fname '_dynamic'],z(iyr0),exo_simul, ...
+                            M.params, dr.ys, it_);
+    end;
+elseif options.order == 2
+    if (options.bytecode)
+        [chck, junk, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ...
+                            M.params, dr.ys, 1);
+        jacobia_ = [loc_dr.g1 loc_dr.g1_x];
+    else
+        [junk,jacobia_,hessian1] = feval([M.fname '_dynamic'],z(iyr0),...
+                                         exo_simul, ...
+                                         M.params, dr.ys, it_);
+    end;
+    if options.use_dll
+        % In USE_DLL mode, the hessian is in the 3-column sparse representation
+        hessian1 = sparse(hessian1(:,1), hessian1(:,2), hessian1(:,3), ...
+                          size(jacobia_, 1), size(jacobia_, 2)*size(jacobia_, 2));
+    end
+end
+
+if any(any(isinf(jacobia_) | isnan(jacobia_)))
+    problem_dummy=1;
+    [infrow,infcol]=find(isinf(jacobia_) | isnan(jacobia_));
+    fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains Inf or NaN. The problem arises from: \n\n')
+    display_problematic_vars_Jacobian(infrow,infcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ')
+end
+if any(any(isinf(hessian1) | isnan(hessian1)))
+    problem_dummy=1;
+    fprintf('\nMODEL_DIAGNOSTICS: The Hessian of the dynamic model contains Inf or NaN.\n')
 end
 
 if problem_dummy==0
-    fprintf('model_diagnostics was not able to detect any obvious problems with this mod-file.\n')
+    fprintf('MODEL_DIAGNOSTICS:  No obvious problems with this mod-file were detected.\n')
 end
 
diff --git a/matlab/print_info.m b/matlab/print_info.m
index 7914043904defb569d193d3bfc1115af17275a40..34017e513f67f1fbb1e23a5824b40b420c03dde1 100644
--- a/matlab/print_info.m
+++ b/matlab/print_info.m
@@ -62,7 +62,7 @@ if ~noprint
           error(['The Jacobian contains NaNs because the following parameters are NaN: '...
               disp_string])
         else
-          error(['The Jacobian contains NaNs'])
+          error(['The Jacobian contains NaNs. For more information, use options_.debug.'])
         end
       case 9
         error(['k_order_pert was unable to compute the solution']) 
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index dd835e2ba611eb90ff732f71e1ccafd0e8406723..2738afc292185c2eee9e08cc99a806c0ead5ee2c 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -119,27 +119,10 @@ end
     
 if options_.debug
     if ~isempty(infrow)     
-        for ii=1:length(infrow)
-            [var_row,var_index]=find(M_.lead_lag_incidence==infcol(ii));
-            if var_row==2
-                type_string='';
-            elseif var_row==1
-                type_string='lag of';
-            elseif var_row==3;
-                type_string='lead of';
-            end
-            if var_index<=M_.orig_endo_nbr
-                fprintf('STOCHASTIC_SOLVER:  Derivative of Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',infrow(ii),type_string,deblank(M_.endo_names(var_index,:)),deblank(M_.endo_names(var_index,:)),dr.ys(var_index))
-            else %auxiliary vars
-                orig_var_index=M.aux_vars(1,var_index-M_.orig_endo_nbr).orig_index;
-                fprintf('STOCHASTIC_SOLVER:  Derivative of Equation %d with respect to %s Variable %s  (initial value of %s: %g) \n',infrow(ii),type_string,deblank(M_.endo_names(orig_var_index,:)),deblank(M_.endo_names(orig_var_index,:)),dr.ys(orig_var_index))            
-            end    
-        end
-        fprintf('\nSTOCHASTIC_SOLVER:  The problem most often occurs, because a variable with\n')
-        fprintf('STOCHASTIC_SOLVER:  exponent smaller than 1 has been initialized to 0. Taking the derivative\n')
-        fprintf('STOCHASTIC_SOLVER:  and evaluating it at the steady state then results in a division by 0.\n')
-    end
+    fprintf('\nSTOCHASTIC_SOLVER: The Jacobian of the dynamic model contains Inf. The problam is associated with:\n\n')    
+    display_problematic_vars_Jacobian(infrow,infcol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ')
     save([M_.fname '_debug.mat'],'jacobia_')
+    end
 end
 
 if ~isempty(infrow)
@@ -157,7 +140,16 @@ if ~isreal(jacobia_)
     end
 end
 
-if any(any(isnan(jacobia_)))
+[nanrow,nancol]=find(isnan(jacobia_));
+if options_.debug
+    if ~isempty(nanrow)     
+    fprintf('\nSTOCHASTIC_SOLVER: The Jacobian of the dynamic model contains NaN. The problam is associated with:\n\n')    
+    display_problematic_vars_Jacobian(nanrow,nancol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ')
+    save([M_.fname '_debug.mat'],'jacobia_')
+    end
+end
+
+if ~isempty(nanrow)
    info(1) = 8;
    NaN_params=find(isnan(M_.params));
    info(2:length(NaN_params)+1) =  NaN_params;
@@ -313,4 +305,5 @@ if options_.loglinear
        error('Loglinear options currently only works at order 1')
     end
 end
+end