diff --git a/equationcollection/equationcollection.py b/equationcollection/equationcollection.py
index c3a6e9294a64abda84692f0207921ba8c9a07643..91e483ef48bee5c4507af449c9e86f796b6e07c4 100644
--- a/equationcollection/equationcollection.py
+++ b/equationcollection/equationcollection.py
@@ -1,6 +1,6 @@
 import sympy as sp
-from copy import copy, deepcopy
-from pystencils.sympyextensions import fastSubs, countNumberOfOperations
+from copy import deepcopy
+from pystencils.sympyextensions import fastSubs, countNumberOfOperations, sortEquationsTopologically
 
 
 class EquationCollection:
@@ -64,7 +64,7 @@ class EquationCollection:
         if addSubstitutionsAsSubexpressions:
             newSubexpressions = [sp.Eq(b, a) for a, b in substitutionDict.items()] + newSubexpressions
 
-        return self.copy(newEquations, newSubexpressions)
+        return self.copy(newEquations, sortEquationsTopologically(newSubexpressions))
 
     def addSimplificationHint(self, key, value):
         """
diff --git a/field.py b/field.py
index 8913fc52a2a35daaff54430d261b46bac6f4dc27..f4424465fbf16f9e46429c9e5507040f59249746 100644
--- a/field.py
+++ b/field.py
@@ -311,7 +311,7 @@ def numpyDataTypeToC(dtype):
         return "float"
     elif dtype == np.int32:
         return "int"
-    raise NotImplementedError()
+    raise NotImplementedError("Cannot convert type " + str(dtype))
 
 
 def offsetComponentToDirectionString(coordinateId, value):
diff --git a/sympyextensions.py b/sympyextensions.py
index a5121775f9c09cd11da3619ba20855e7d7503cf0..1d3203e662572e8a30181385d529d1ea9193c30c 100644
--- a/sympyextensions.py
+++ b/sympyextensions.py
@@ -367,3 +367,8 @@ def getSymmetricPart(term, vars):
     """
     substitutionDict = {e: -e for e in vars}
     return sp.Rational(1, 2) * (term + fastSubs(term, substitutionDict))
+
+
+def sortEquationsTopologically(equationSequence):
+    res = sp.cse_main.reps_toposort([[e.lhs, e.rhs] for e in equationSequence])
+    return [sp.Eq(a, b) for a, b in res]