diff --git a/pymule/colours.py b/pymule/colours.py
new file mode 100644
index 0000000000000000000000000000000000000000..5247cb7614a52978f04d3a88c080ece99ce7b3ad
--- /dev/null
+++ b/pymule/colours.py
@@ -0,0 +1,25 @@
+from cycler import cycler
+from matplotlib import rc
+
+
+# Mathematica colour scheme 112 `vibrant'
+
+# Apply[
+#   StringJoin["#", ##] &,
+#   Map[
+#     IntegerString[#, 16, 2] &,
+#     Round[255 List @@@ ColorData[112] /@ {2, 3, 4, 1, 5, 6, 7, 8, 9, 10}]
+#   ],
+#   {1}
+# ]
+schema = [
+  "#3163ce", "#ff9b00", "#00981c", "#ca3300", "#9152ba",
+  "#0096b4", "#d96d22", "#7f49c6", "#20a178", "#d24b31"
+]
+rc('axes', prop_cycle=cycler(u'color', schema))
+
+orderscheme = {
+  'lo': 'C2', 'nlo': 'C0', 'nnlo': 'C1', 'np': 'C3'
+}
+
+defcol = ['C%d' % i for i in range(10)]
diff --git a/pymule/plot.py b/pymule/plot.py
index d07918e98fdfb0dc47a6ee65b12cf4bdb1dcceff..77bc629ca553f896ab1ec52d26dcf6cbed95cef7 100644
--- a/pymule/plot.py
+++ b/pymule/plot.py
@@ -3,6 +3,7 @@ import matplotlib.collections
 import matplotlib.patches
 import numpy as np
 from matplotlib import rc
+from colours import *
 from errortools import *
 
 
@@ -13,7 +14,7 @@ rc('text.latex', preamble="\n".join([
 ]))
 
 
-def errorband(p, ax=None, underflow=False, overflow=False):
+def errorband(p, ax=None, col='default', underflow=False, overflow=False):
     if ax is None:
         ax = plt.gca()
 
@@ -41,8 +42,11 @@ def errorband(p, ax=None, underflow=False, overflow=False):
         for x, y, e in p
     ]
 
-    artist = ax.step(p[:,0], p[:,1])
-    col = artist[0].get_color()
+    if col == 'default':
+        artist = ax.step(p[:,0], p[:,1])
+        col = artist[0].get_color()
+    else:
+        artist = ax.step(p[:,0], p[:,1], col)
 
     pc = matplotlib.collections.PatchCollection(
         boxes,
@@ -54,10 +58,10 @@ def errorband(p, ax=None, underflow=False, overflow=False):
 
 
 def twopanel(labelx,
-             upleft=[], labupleft="",
-             downleft=[], labdownleft="",
-             upright=[], labupright="",
-             downright=[], labdownright=""):
+             upleft=[], labupleft="", colupleft=defcol,
+             downleft=[], labdownleft="", coldownleft=defcol,
+             upright=[], labupright="", colupright=defcol,
+             downright=[], labdownright="", coldownright=defcol):
 
     if type(upleft) == np.ndarray:
         upleft = [upleft]
@@ -74,26 +78,26 @@ def twopanel(labelx,
     axs[1].set_xlabel(labelx)
 
     axs[0].set_ylabel(labupleft)
-    for i in upleft:
-        errorband(i, ax=axs[0])
+    for i, c in zip(upleft, colupleft):
+        errorband(i, ax=axs[0], col=c)
 
     if len(upright) > 0:
         ax2 = axs[0].twinx()
         if labupright is not None:
             ax2.set_ylabel(labupright)
-        for i in upright:
-            errorband(i, ax=ax2)
+        for i, c in zip(upright, colupright):
+            errorband(i, ax=ax2, col=c)
 
     axs[1].set_ylabel(labdownleft)
-    for i in downleft:
-        errorband(i, ax=axs[1])
+    for i,c  in zip(downleft, coldownleft):
+        errorband(i, ax=axs[1], col=c)
 
     if len(downright) > 0:
         ax3 = axs[1].twinx()
         if labdownright is not None:
             ax3.set_ylabel(labdownright)
-        for i, c in downright:
-            errorband(i, ax=ax3)
+        for i, c in zip(downright, coldownright):
+            errorband(i, ax=ax3, col=c)
 
     return fig, axs
 
@@ -135,11 +139,16 @@ def kplot(sigma, labelx='x_e', labelsigma=None,
 
     if 1 in showk:
         kwargs['downleft'] = divideplots(sigma['nlo'], sigma['lo'])
+        kwargs['coldownleft'] = [orderscheme['nlo']]
     if 2 in showk and 'nnlo' in sigma:
         kwargs['downright'] = divideplots(sigma['nnlo'], sigma['lo'])
+        kwargs['coldownright'] = [orderscheme['nnlo']]
 
     kwargs['upleft'] = [
         xsec[orders[i]] for i in show
     ]
+    kwargs['colupleft'] = [
+        orderscheme[orders[i]] for i in show
+    ]
 
     return twopanel(**kwargs)