diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 76d1be395373a7c6c82e1f24ea4b4c4e60c69a11..c59b9be10e68b3817a0c507eb3f3d450a6293d0e 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -4194,6 +4194,19 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
   // Computes dynamic jacobian columns, must be done after computeDerivIDs()
   computeDynJacobianCols(jacobianExo);
 
+  /* In both MATLAB and Julia, tensors for higher-order derivatives are stored
+     in matrices whose columns correspond to variable multi-indices. Since we
+     currently are limited to 32-bit signed integers (hence 31 bits) for matrix
+     indices, check that we will not overflow (see #89). Note that such a check
+     is not needed for parameter derivatives, since tensors for those are not
+     stored as matrices. This check cannot be done before since
+     dynJacobianColsNbr is not yet set.*/
+  if (log2(dynJacobianColsNbr)*derivsOrder >= numeric_limits<int>::digits)
+    {
+      cerr << "ERROR: The dynamic derivatives matrix is too large. Please decrease the approximation order." << endl;
+      exit(EXIT_FAILURE);
+    }
+
   // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic
   set<int> vars;
   for (auto &it : deriv_id_table)
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 33973cf95bffa6c597166cd2c46bb3b6247e7e5d..906bdca0d7ec919f152bd8be5cadf16940f83214 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1003,6 +1003,19 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
   equations.clear();
   copy(neweqs.begin(), neweqs.end(), back_inserter(equations));
 
+  /* In both MATLAB and Julia, tensors for higher-order derivatives are stored
+     in matrices whose columns correspond to variable multi-indices. Since we
+     currently are limited to 32-bit signed integers (hence 31 bits) for matrix
+     indices, check that we will not overflow (see #89). Note that such a check
+     is not needed for parameter derivatives, since tensors for those are not
+     stored as matrices. This check is implemented at this place for symmetry
+     with DynamicModel::computingPass(). */
+  if (log2(symbol_table.endo_nbr())*derivsOrder >= numeric_limits<int>::digits)
+    {
+      cerr << "ERROR: The static derivatives matrix is too large. Please decrease the approximation order." << endl;
+      exit(EXIT_FAILURE);
+    }
+
   // Compute derivatives w.r. to all endogenous
   set<int> vars;
   for (int i = 0; i < symbol_table.endo_nbr(); i++)