diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 5db9b7bcb46aedd4d95c252f74508ee3bd982e8c..2c2f41dfd2279470f4e3b20ec944e7f6c619c25a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -19,4 +19,4 @@ target_compile_definitions("mmdg-3d" PRIVATE GRIDDIM=3)
 target_link_dune_default_libraries("mmdg-3d")
 
 dune_symlink_to_source_files(FILES grids parameterMMDG.ini parameterDG.ini
-  dgAnalysis.py mmdgAnalysis.py)
+  dgAnalysis.py mmdgAnalysis.py plotDG.py plotMMDG.py)
diff --git a/src/plotDG.py b/src/plotDG.py
new file mode 100644
index 0000000000000000000000000000000000000000..dbbed59eac9b5503802f0b16e6d93c287bfe3bbb
--- /dev/null
+++ b/src/plotDG.py
@@ -0,0 +1,30 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpltools.annotation import slope_marker
+import matplotlib.pylab as pylab
+
+params = {'legend.fontsize': 'x-large',
+         'axes.labelsize': 'x-large',
+         'axes.titlesize':'x-large',
+         'xtick.labelsize':'x-large',
+         'ytick.labelsize':'x-large'}
+pylab.rcParams.update(params)
+
+data_1d = np.loadtxt("plots/dgAnalysis_1d")
+data_2d = np.loadtxt("plots/dgAnalysis_2d")
+data_3d = np.loadtxt("plots/dgAnalysis_3d")
+
+plt.loglog(np.reciprocal(data_1d[:,2]), data_1d[:,0], 'g^-', label='$n=1$')
+plt.loglog(np.reciprocal(data_2d[:,2]), data_2d[:,0], 'bo-', label='$n=2$')
+plt.loglog(np.reciprocal(data_3d[:,2]), data_3d[:,0], 'rs-', label='$n=3$')
+
+x = 1.0 / data_3d[3,2]
+y = 10 ** (0.2 * np.log10(data_3d[3,0]) + 0.8 * np.log10(data_3d[4,0]))
+slope_marker((x, y), -2, size_frac = 0.15, invert = True, pad_frac = 0.12, \
+    poly_kwargs={'ec': 'black', 'fill': False}, \
+    text_kwargs={'size': 'x-large', 'text': '2', 'usetex': True})
+plt.xlabel('$h^{-1}$')
+plt.ylabel('$\Vert p - p_h \Vert_{L^2 ( \Omega )}$')
+plt.legend()
+plt.tight_layout()
+plt.savefig('plots/dg.pdf', format='pdf')
diff --git a/src/plotMMDG.py b/src/plotMMDG.py
new file mode 100644
index 0000000000000000000000000000000000000000..40a6e68899afebaf702e069f05d59e7dd0b0dd9a
--- /dev/null
+++ b/src/plotMMDG.py
@@ -0,0 +1,68 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpltools.annotation import slope_marker
+import matplotlib.pylab as pylab
+
+params = {'legend.fontsize': 'x-large',
+         'axes.labelsize': 'x-large',
+         'axes.titlesize':'x-large',
+         'xtick.labelsize':'x-large',
+         'ytick.labelsize':'x-large'}
+pylab.rcParams.update(params)
+
+data_2d = np.loadtxt("plots/mmdgAnalysis_2d")
+data_3d = np.loadtxt("plots/mmdgAnalysis_3d")
+
+# plot total error
+plt.figure()
+plt.loglog(np.reciprocal(np.maximum(data_2d[:,3], data_2d[:,4])), data_2d[:,2],\
+    'bo-', label='$n=2$')
+plt.loglog(np.reciprocal(np.maximum(data_3d[:,3], data_3d[:,4])), data_3d[:,2],\
+    'rs-', label='$n=3$')
+
+if data_3d[0,2] > data_3d[-1,2]:
+    x = 1.0 / np.maximum(data_3d[3,3], data_3d[3,4])
+    y = 10 ** (0.2 * np.log10(data_3d[3,2]) + 0.8 * np.log10(data_3d[4,2]))
+    slope_marker((x,y), -2, size_frac = 0.15, invert = True, pad_frac = 0.12, \
+        poly_kwargs={'ec': 'black', 'fill': False}, \
+        text_kwargs={'size': 'x-large', 'text': '2', 'usetex': True})
+plt.xlabel('$h^{-1}$')
+plt.ylabel(r'$\Vert (p,p^\Gamma) - (p_h , p_h^\Gamma) \Vert_{L^2 ( \Omega ) '\
+    + r'\times L^2 ( \Gamma )}$')
+plt.legend()
+plt.tight_layout()
+plt.savefig('plots/mmdg_total.pdf', format='pdf')
+
+# plot bulk error
+plt.figure()
+plt.loglog(np.reciprocal(data_2d[:,3]), data_2d[:,0], 'bo-', label='$n=2$')
+plt.loglog(np.reciprocal(data_3d[:,3]), data_3d[:,0], 'rs-', label='$n=3$')
+
+if data_3d[0,0] > data_3d[-1,0]:
+    x = 1.0 / data_3d[3,3]
+    y = 10 ** (0.2 * np.log10(data_3d[3,0]) + 0.8 * np.log10(data_3d[4,0]))
+    slope_marker((x,y), -2, size_frac = 0.15, invert = True, pad_frac = 0.12, \
+        poly_kwargs={'ec': 'black', 'fill': False}, \
+        text_kwargs={'size': 'x-large', 'text': '2', 'usetex': True})
+plt.xlabel('$h^{-1}$')
+plt.ylabel(r'$\Vert p - p_h  \Vert_{L^2 ( \Omega )}$')
+plt.legend()
+plt.tight_layout()
+plt.savefig('plots/mmdg_bulk.pdf', format='pdf')
+
+# plot interface error
+plt.figure()
+plt.loglog(np.reciprocal(data_2d[:,4]), data_2d[:,1], 'bo-', label='$n=2$')
+plt.loglog(np.reciprocal(data_3d[:,4]), data_3d[:,1], 'rs-', label='$n=3$')
+
+if data_3d[0,1] > data_3d[-1,1]:
+    x = 1.0 / data_3d[3,4]
+    y = 10 ** (0.1 * np.log10(data_3d[3,1]) + 0.9 * np.log10(data_3d[4,1]))
+    slope_marker((x,y), -2, size_frac = 0.15, invert = True, pad_frac = 0.12, \
+        poly_kwargs={'ec': 'black', 'fill': False}, \
+        text_kwargs={'size': 'x-large', 'text': '2', 'usetex': True})
+plt.xlabel('$h^{-1}$')
+plt.ylabel(r'$\Vert p^\Gamma - p_h^\Gamma \Vert_{L^2 ( \Gamma )}$')
+plt.legend()
+plt.tight_layout()
+plt.savefig('plots/mmdg_interface.pdf', format='pdf')