diff --git a/jastadd-mquat-benchmark/results/plot.py b/jastadd-mquat-benchmark/results/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..e81234f646185f92c82315e29bfed447d7faef0d
--- /dev/null
+++ b/jastadd-mquat-benchmark/results/plot.py
@@ -0,0 +1,153 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from matplotlib import colors as mcolors
+
+# constants
+bar_width = 0.2
+solver_names = ['ilp-direct', 'ilp-external', 'random', 'simple']
+patterns = ['❌', '✔']
+
+colors = [c[4:] for c in (sorted(mcolors.TABLEAU_COLORS.keys())) if 'dark' + c[4:] in mcolors.CSS4_COLORS]
+colors.reverse()
+
+
+def to_label(name):
+    if '-' in name:
+        tokens = name.split('-')
+        return '{} ({})'.format(tokens[0].upper(), tokens[1].title())
+    else:
+        return name.title()
+
+
+def load(name, show_head=False):
+    data = pd.read_csv(name)
+    data['tool'] = data.name.astype(str)
+    data['Gen'] = (data.Gen.astype(int) > 0) * data.Gen.astype(int)
+    data['total'] = data.Gen.astype(int) + data.Solved.astype(int)
+    m = data.groupby(['name'])[['total']].apply(np.median)
+    m.name = 'MEDIAN'
+    data = data.join(m, on=['name'])
+    data = data.drop_duplicates(subset='tool')
+    data = data[data.tool.isin(solver_names)]
+    if show_head:
+        data.head()
+    return data
+
+
+def load_scenario_file(name, show_head=False):
+    for filename in os.listdir("scenarios"):
+        if filename.endswith(name):
+            return load("scenarios/" + filename, show_head)
+
+
+def plot_range(names, labels, ax):
+
+    # constants
+    N = len(solver_names)
+    width = 0.05
+    ind = np.arange(N)
+
+    labelpos = []
+    for i in range(len(names)):
+        labelpos += [(N/2 - 0.5 + i * (N + 0.5)) * width]
+
+    for (name, i) in zip(names, range(len(names))):
+
+        # load data
+        data = load_scenario_file(name)
+
+        print(data)
+
+
+        means = data.MEDIAN
+
+
+        rect = ax.bar((ind + i * (N + 0.5)) * width, means, width=width, color=colors)
+        for r, isValid, color in zip(rect, data.Valid, colors):
+            ax.text(r.get_x() + r.get_width() / 2,
+                    r.get_height() * 1.5,
+                    patterns[isValid],  # the text
+                    fontname='symbola', fontsize=16, ha='center', va='bottom', color='black')
+
+        for r, isValid, color in zip(rect, data.Valid, colors):
+            ax.text(r.get_x() + r.get_width() / 2,
+                    r.get_height(),
+                    '%d' % int(r.get_height()),
+                    fontsize=8, ha='center', va='bottom', color='black')
+
+    ax.set_xticks(labelpos)
+    ax.set_xticklabels(labels)
+    ax.set_ylim([1, 6000000])
+
+    ax.set_yticks([1, 1000, 1000000])
+
+    return rect
+
+
+def create_grid():
+    variants = [2, 4]
+    requests = [1, 2, 3]
+    depth = [1, 2, 3]
+    resources = [15, 30]
+
+    y_dimension = (requests, "Requests")
+    x_dimension = (depth, "Depth")
+    inner_dimension = (variants, "Variants")
+    fixed_dimension = (resources, "Resources")
+
+    pos_map = {
+        "Variants": 0,
+        "Requests": 1,
+        "Depth": 2,
+        "Resources": 3
+    }
+
+    for fixed_pos in fixed_dimension[0]:
+        fig, axs = plt.subplots(len(y_dimension[0]), len(x_dimension[0]), figsize=(10, 10), subplot_kw=dict(yscale="log"), gridspec_kw={})
+        for y, y_pos in zip(y_dimension[0], range(len(y_dimension[0]))):
+            for x, x_pos in zip(x_dimension[0], range(len(x_dimension[0]))):
+                parameters = []
+                labels = []
+                for inner in inner_dimension[0]:
+                    new_parameters = [0,0,0,0]
+                    new_parameters[pos_map[x_dimension[1]]] = x
+                    new_parameters[pos_map[y_dimension[1]]] = y
+                    new_parameters[pos_map[inner_dimension[1]]] = inner
+                    new_parameters[pos_map[fixed_dimension[1]]] = fixed_pos
+                    parameters += [new_parameters]
+                    labels += [inner_dimension[1] + ' = %d' % inner]
+                rect = plot_range(['size_v%d_q%d_d%d_r%d.csv' % (p[0], p[1], p[2], p[3]) for p in parameters], [label for label in labels], axs[y_pos, x_pos])
+                if x_pos == 0 and y_pos == 0:
+                    axs[y_pos, x_pos].set_ylabel("X",fontsize=70,color="white")
+
+                if x_pos == 0 and y_pos ==  len(y_dimension[0]) - 1:
+                    axs[y_pos, x_pos].set_xlabel("X",fontsize=70,color="white")
+                    axs[y_pos, x_pos].legend(rect, [to_label(n) for n in solver_names],
+                                             bbox_to_anchor=(0., -0.3, len(x_dimension[0]), 0.), loc=0, ncol=2, borderaxespad=0.)
+
+        # set the title of the individual rows
+        for i in range(len(y_dimension[0])):
+            axs[i][-1].set_ylabel(y_dimension[1] + ' = ' + str(y_dimension[0][i]), rotation=90)
+            axs[i][-1].yaxis.set_label_coords(1.2, 0.5)
+
+        # set the titles of the individual columns (and enlarge the canvas to make space for the title)
+        for i in range(len(x_dimension[0])):
+            axs[0][i].set_xlabel("\n\n\n\n" + x_dimension[1] + ' = ' + str(x_dimension[0][i]), va="center")
+            axs[0][i].xaxis.set_label_coords(0.5, 1.2)
+
+        # set the background color to white (only required in pycharm)
+        fig.patch.set_facecolor('w')
+
+        # label the common y axis
+        fig.text(0.06, 0.5, 'Runtime [ms]', ha='center', va='center', rotation='vertical')
+
+        # set the figure title to the
+        fig.suptitle("%s = %d" % (fixed_dimension[1], fixed_pos), fontsize=14)
+
+        plt.show()
+
+
+create_grid()