diff --git a/exercises/exercise-coupling-ff-pm/models/plotFluxes.py b/exercises/exercise-coupling-ff-pm/models/plotFluxes.py
new file mode 100644
index 0000000000000000000000000000000000000000..832062f0fc2b4c9945f22d1c384cb413a2f7d151
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/plotFluxes.py
@@ -0,0 +1,32 @@
+import json
+import matplotlib.pyplot as plt
+
+# JSON file
+jsonFilePath = "./flux_models_coupling.json"
+file = open(jsonFilePath)
+fullEvaporationData = json.loads(file.read())
+
+for key in fullEvaporationData:
+    print("developing plot " + key)
+
+    fig, ax = plt.subplots(1, 1, figsize=(6, 6))
+    plt.title('Interfacial Evaporation Rate: ' + key, fontsize=24)
+    plt.xlabel('Interface x')
+    plt.ylabel('Evaporation flux')
+    evaporationData = fullEvaporationData[key]
+    coordinates = evaporationData["FaceCenterCoordinates"]
+    timeSeries = evaporationData["Time"]
+    xCoords = [x for x,y in coordinates]
+    facewiseEvaporationResults = evaporationData["FaceEvaporation"]
+
+    count = 0
+    for time, results in zip(timeSeries, facewiseEvaporationResults):
+        if (count%2 == 0 and time >= 1.0): # plot every 5th time step
+            ax.plot(xCoords, results, label="{:0.3f}".format(time))
+        count = count+1
+
+    ax.legend(title="Time (s)")
+    outputFig = plt.gcf()
+    outputFig.savefig('./EvaporationOverInterface_' + key + '.png', bbox_inches='tight')
+    plt.clf()
+    plt.close()
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py b/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py
new file mode 120000
index 0000000000000000000000000000000000000000..a3955ceebbe41d3cc49b61f87d939524021cf5bd
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/plotFluxes.py
@@ -0,0 +1 @@
+/home/nedc/DUMUX_summerschool/dumux-course/exercises/exercise-coupling-ff-pm/models/plotFluxes.py
\ No newline at end of file