diff --git a/dagrt/language.py b/dagrt/language.py
index a9f5dc1a5c4c494ac295828ad3b8df255559d36f..60a3709a694ba58924945d0ca7c362c029f33f5a 100644
--- a/dagrt/language.py
+++ b/dagrt/language.py
@@ -31,7 +31,7 @@ from pymbolic.imperative.statement import (
         ConditionalAssignment as AssignBase,
         Nop as NopBase)
 
-from dagrt.utils import get_variables
+from dagrt.utils import get_variables, natsorted
 from contextlib import contextmanager
 
 import logging
@@ -162,7 +162,7 @@ def _stringify_statements(roots, id_to_stmt, prefix=""):
             return
         printed_stmt_ids.add(stmt.id)
 
-        for dep_id in stmt.depends_on:
+        for dep_id in natsorted(stmt.depends_on):
             print_stmt(id_to_stmt[dep_id])
 
         lines.append(
@@ -693,14 +693,17 @@ class DAGCode(RecordWithoutPickling):
     def __str__(self):
         lines = []
         for phase_name, phase in sorted(six.iteritems(self.phases)):
-            lines.append("STAGE %s" % phase_name)
+            phase_title = "PHASE \"%s\"" % phase_name
+            if phase_name == self.initial_phase:
+                phase_title += " (initial_phase)"
+            lines.append(phase_title)
 
-            for root_id in phase.depends_on:
+            for root_id in natsorted(phase.depends_on):
                 lines.extend(_stringify_statements(
                     [phase.id_to_stmt[root_id]],
                     phase.id_to_stmt, prefix="    "))
 
-            lines.append("    -> (next phase) %s" % phase.next_phase)
+            lines.append("    -> (next phase) \"%s\"" % phase.next_phase)
             lines.append("")
 
         return "\n".join(lines)
@@ -1165,14 +1168,14 @@ def get_dot_dependency_graph(code, use_stmt_ids=False):
     def additional_lines_hook():
         for i, (name, phase) in enumerate(six.iteritems(code.phases)):
             yield "subgraph cluster_%d { label=\"%s\"" % (i, name)
-            for dep in phase.depends_on:
+            for dep in natsorted(phase.depends_on):
                 yield dep
             yield "}"
 
     statements = [
             stmt if use_stmt_ids else stmt.copy(id=stmt.id)
             for phase_name, phase in six.iteritems(code.phases)
-            for stmt in phase.statements]
+            for stmt in natsorted(phase.statements, key=lambda stmt: stmt.id)]
 
     return get_dot_dependency_graph(
             statements, use_stmt_ids=use_stmt_ids,
diff --git a/dagrt/utils.py b/dagrt/utils.py
index bd8df2e9fbec232b721e49eb4d59d66a4c0b2500..69f4640936b1abeba6e66851e1c2eb5d8ab6c0b6 100644
--- a/dagrt/utils.py
+++ b/dagrt/utils.py
@@ -224,4 +224,33 @@ def run_fortran(sources, fortran_options=None, fortran_libraries=None):
 
 # }}}
 
+
+# {{{ sorting in natural order
+
+def natorder(key):
+    # Return natural ordering for strings, as opposed to dictionary order.
+    # E.g. will result in
+    #  'abc1' < 'abc9' < 'abc10'
+    # rather than
+    #  'abc1' < 'abc10' < 'abc9'
+    # Based on
+    # http://code.activestate.com/recipes/285264-natural-string-sorting/#c7
+    import re
+    result = []
+    for (int_val, string_val) in re.findall(r"(\d+)|(\D+)", key):
+        if int_val:
+            result.append(int(int_val))
+            # Tie-breaker in case leading zeros in *int_val* cause distinct
+            # values to compare equally.
+            result.append(len(int_val))
+        else:
+            result.append(string_val)
+    return result
+
+
+def natsorted(seq, key=lambda x: x):
+    return sorted(seq, key=lambda y: natorder(key(y)))
+
+# }}}
+
 # vim: foldmethod=marker