diff --git a/matlab/print_info.m b/matlab/print_info.m
index a4b775d39ebc2dc101514c891f3657b4a1ee0a54..0ec5e378884827ee3e9790ad901c4fb9c1db2e3f 100644
--- a/matlab/print_info.m
+++ b/matlab/print_info.m
@@ -68,6 +68,10 @@ if ~noprint
         error(['k_order_pert was unable to compute the solution'])
       case 10
         error(['The Jacobian of the dynamic model contains Inf. For more information, use options_.debug.'])
+      case 11
+        error('The Hessian of the dynamic model used for second order solutions must not contain Inf')
+      case 12
+        error('The Hessian of the dynamic model used for second order solutions must not contain NaN')
       case 19
         error('The steadystate file did not compute the steady state')
       case 20
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index 2738afc292185c2eee9e08cc99a806c0ead5ee2c..8f22a025ce62717203399927590206bba79dfc41 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -113,6 +113,30 @@ elseif options_.order == 2
         hessian1 = sparse(hessian1(:,1), hessian1(:,2), hessian1(:,3), ...
                           size(jacobia_, 1), size(jacobia_, 2)*size(jacobia_, 2));
     end
+    [infrow,infcol]=find(isinf(hessian1));
+    if options_.debug
+        if ~isempty(infrow)     
+        fprintf('\nSTOCHASTIC_SOLVER: The Hessian of the dynamic model contains Inf.\n')
+        fprintf('STOCHASTIC_SOLVER: Try running model_diagnostics to find the source of the problem.\n')
+        save([M_.fname '_debug.mat'],'hessian1')
+        end
+    end
+    if ~isempty(infrow)
+        info(1)=11;
+        return
+    end
+    [nanrow,nancol]=find(isnan(hessian1));
+    if options_.debug
+        if ~isempty(nanrow)     
+            fprintf('\nSTOCHASTIC_SOLVER: The Hessian of the dynamic model contains NaN.\n')
+            fprintf('STOCHASTIC_SOLVER: Try running model_diagnostics to find the source of the problem.\n')
+            save([M_.fname '_debug.mat'],'hessian1')
+        end
+    end
+    if ~isempty(nanrow)
+        info(1)=12;
+        return
+    end    
 end
 
 [infrow,infcol]=find(isinf(jacobia_));