diff --git a/harmonicmaps-skyrmions-eoc-testsuite/testsuite.py b/harmonicmaps-skyrmions-eoc-testsuite/testsuite.py
new file mode 100644
index 0000000000000000000000000000000000000000..0b186e6e1dfb85fe5fedc03cb05e9099c0fbbf55
--- /dev/null
+++ b/harmonicmaps-skyrmions-eoc-testsuite/testsuite.py
@@ -0,0 +1,36 @@
+import subprocess
+
+for order in range(1,4):
+
+    minLevel =  5 + 1 - order;
+    maxLevel = 10 + 1 - order;
+
+    processList = []
+
+    for numLevels in range(minLevel,maxLevel+1):
+        subprocess.call(["echo", "foo"])
+
+        LOGFILE = "./harmonicmaps_" + str(order) + "_" + str(numLevels) + ".log"
+
+        ##  run the actual simulation
+        executable = "./harmonicmaps-" + str(order)
+        p = subprocess.Popen([executable, "harmonicmaps-skyrmions-hexagon.parset", "-numLevels", str(numLevels)])
+        processList.append(p)
+
+    exit_codes = [p.wait() for p in processList]
+
+    subprocess.call(["echo", "Now measuring errors"])
+
+
+    for numLevels in range(minLevel,maxLevel+1):
+        # Measure the discretization errors against the solution on the finest grid
+        LOGFILE = "./compute-disc-error_" + str(order) + "_" + str(numLevels) + ".log"
+
+        subprocess.Popen(["../build-cmake/src/compute-disc-error", "compute-disc-error-skyrmions-hexagon.parset",
+                                        "-order", str(order),
+                                        "-numLevels", str(numLevels),
+                                        "-numReferenceLevels", str(maxLevel),
+                                        "-simulationData", "harmonicmaps-result-" + str(order) + "-" + str(numLevels) + ".data",
+                                        "-referenceData", "harmonicmaps-result-" + str(order) + "-" + str(maxLevel) + ".data"])
+
+