diff --git a/benchmarks/benchmarks/scripts.py b/benchmarks/benchmarks/scripts.py
index 5202dfd..06bb436 100644
--- a/benchmarks/benchmarks/scripts.py
+++ b/benchmarks/benchmarks/scripts.py
@@ -15,7 +15,7 @@
 # def circuit_benchmark(nqubits, backend, circuit_name, circuit_options=None,
 def circuit_benchmark(
     n_hydrogens,
-    backend,  # circuit_name, circuit_options=None,
+    backend,
     nreps=1,
     nshots=None,
     transfer=False,
diff --git a/benchmarks/data/benchmark.png b/benchmarks/data/benchmark.png
new file mode 100644
index 0000000..9258d31
Binary files /dev/null and b/benchmarks/data/benchmark.png differ
diff --git a/benchmarks/data/plot-graphs.py b/benchmarks/data/plot-graphs.py
new file mode 100644
index 0000000..fbc8df6
--- /dev/null
+++ b/benchmarks/data/plot-graphs.py
@@ -0,0 +1,50 @@
+"""
+Parse the .dat files from running main.py to get and plot the dry_run_times
+"""
+
+import json
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def parse_dat_files(dat_glob, directory=None):
+    """
+    Read all the .dat files found with dat_glob. .dat files should be in .json format
+    """
+    if directory is None:
+        directory = "."
+
+    result = {}
+    for dat_file in Path(directory).glob(dat_glob):
+        _h_n, backend = dat_file.stem.split("_")
+        h_n = int(_h_n[1:])
+        if result.get(h_n) is None:
+            result[h_n] = {}
+        with open(dat_file) as file_handler:
+            result[h_n][backend] = json.load(file_handler)
+    return result
+
+
+def main():
+    dat_glob = "*.dat"
+    json_data = parse_dat_files(dat_glob)
+    # Plot results
+    x_vals = sorted(json_data.keys())
+    # print(x_vals)
+    qibojit_vals = [json_data[_i]["qibojit"][0]["dry_run_time"] for _i in x_vals]
+    # print(qibojit_vals)
+    numpy_vals = [json_data[_i]["numpy"][0]["dry_run_time"] for _i in x_vals]
+    # print(numpy_vals)
+
+    plt.plot(x_vals, qibojit_vals, label="Qibojit")
+    plt.plot(x_vals, numpy_vals, label="Numpy")
+    plt.legend()
+    plt.ylabel("Time (s)")
+    plt.xlabel(r"Length of hydrogen chain, $H_n$")
+    plt.savefig("benchmark.svg")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/benchmarks/run-benchmark.sh b/benchmarks/run-benchmark.sh
new file mode 100755
index 0000000..5122207
--- /dev/null
+++ b/benchmarks/run-benchmark.sh
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+
+nreps=10
+
+# Run main.py with various command-line arguments/options
+for n_hydrogen in 2 4 6 8; do
+    for backend in numpy qibojit; do
+        filename="h${n_hydrogen}_${backend}"
+        for precision in single double; do
+            python main.py --n_hydrogens ${n_hydrogen} --backend ${backend} --nreps ${nreps} --precision ${precision} --filename ${filename}.dat &>> ${filename}.log
+            echo &>> ${filename}.log
+        done
+
+        # Move the .log and .dat file into data/ folder
+        if [ -d "data" ]; then
+            mv ${filename}.dat ${filename}.log
+        fi
+    done
+done