diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 6e8fcdfebf0250cbb65c3720c414dc40aee1f6cc..03a3b020b2f6d4bc238cd8bdbef0aad6107ae23e 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -1741,7 +1741,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
         {
           // Writing initialization instruction for matrix g2
           DynamicOutput << "if nargout >= 3," << endl
-                        << "  g2 = sparse([],[],[], " << nrows << ", " << hessianColsNbr << ", " << 5*hessianColsNbr << ");" << endl
+                        << "  g2 = sparse([],[],[], " << nrows << ", " << hessianColsNbr << ", " << NNZDerivatives[1] << ");" << endl
                         << endl
                         << "%" << endl
                         << "% Hessian matrix" << endl
@@ -1755,7 +1755,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
         {
           int ncols = hessianColsNbr * dynJacobianColsNbr;
           DynamicOutput << "if nargout >= 4," << endl
-                        << "  g3 = sparse([],[],[], " << nrows << ", " << ncols << ", " << 5*ncols << ");" << endl
+                        << "  g3 = sparse([],[],[], " << nrows << ", " << ncols << ", " << NNZDerivatives[2] << ");" << endl
                         << endl
                         << "%" << endl
                         << "% Third order derivatives" << endl
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index ed7f578c1f91504cf5e374517923661d26b5de6a..ddb302a4963cb5baad88d0efa2e610cfe5f0df72 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -28,6 +28,8 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg,
   DataTree(symbol_table_arg, num_constants_arg),
   mode(eStandardMode)
 {
+  for(int i=0; i < 3; i++)
+    NNZDerivatives[i] = 0;
 }
 
 int
@@ -59,6 +61,7 @@ ModelTree::computeJacobian(const set<int> &vars)
         if (d1 == Zero)
           continue;
         first_derivatives[make_pair(eq, *it)] = d1;
+	++NNZDerivatives[0];
       }
 }
 
@@ -84,6 +87,10 @@ ModelTree::computeHessian(const set<int> &vars)
           if (d2 == Zero)
             continue;
           second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
+	  if (var2 == var1)
+	    ++NNZDerivatives[1];
+	  else
+	    NNZDerivatives[1] += 2;
         }
     }
 }
@@ -114,6 +121,12 @@ ModelTree::computeThirdDerivatives(const set<int> &vars)
           if (d3 == Zero)
             continue;
           third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
+	  if (var3 == var2 && var2 == var1)
+	    ++NNZDerivatives[2];
+	  else if (var3 == var2 || var2 == var1)
+	    NNZDerivatives[2] += 3;
+	  else
+	    NNZDerivatives[2] += 6;
         }
     }
 }
diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh
index 8c61f21d98fb165e3a466f0c91fa012a5839dc85..3d82b4693c3817555e8e38b0989d12d71d155192 100644
--- a/preprocessor/ModelTree.hh
+++ b/preprocessor/ModelTree.hh
@@ -45,6 +45,9 @@ protected:
   //! Stores declared equations
   vector<BinaryOpNode *> equations;
 
+  //! Number of non-zero derivatives
+  int NNZDerivatives[3];
+
   typedef map<pair<int, int>, NodeID> first_derivatives_type;
   //! First order derivatives
   /*! First index is equation number, second is variable w.r. to which is computed the derivative.