diff --git a/exercises/exercise-fluidsystem/2pproblem.hh b/exercises/exercise-fluidsystem/2pproblem.hh
index dc7347aec91ef1b3c109034ef085357e1fa2f89e..0e5d7438aabeb190c60f5721457fd664b38a50a7 100644
--- a/exercises/exercise-fluidsystem/2pproblem.hh
+++ b/exercises/exercise-fluidsystem/2pproblem.hh
@@ -59,6 +59,9 @@
 // The two-phase immiscible fluid system
 #include <dumux/material/fluidsystems/2pimmiscible.hh>
 
+// The interface to create plots during simulation
+#include <dumux/io/gnuplotinterface.hh>
+
 namespace Dumux{
 // Forward declaration of the problem class
 template <class TypeTag> class ExerciseFluidsystemProblemTwoP;
@@ -156,6 +159,8 @@ public:
 
         // set the depth of the bottom of the reservoir
         depthBOR_ = this->fvGridGeometry().bBoxMax()[dimWorld-1];
+        if(getParam<bool>("Output.PlotDensity"))
+            plotDensity_();
     }
 
     /*!
@@ -282,8 +287,27 @@ public:
     }
 
 private:
+    void plotDensity_()
+    {
+        std::vector<double> x(100);
+        std::vector<double> y(100);
+        for (int i=0; i<100; ++i)
+            x[i] = 1e4 + i*100910;
+        for (int i=0; i<100; ++i)
+            y[i] = MyCompressibleComponent<Scalar>::liquidDensity(300, x[i]);
+        gnuplot_.resetPlot();
+        gnuplot_.setXRange(1e4, 1e7);
+        gnuplot_.setOption("set logscale x 10");
+        gnuplot_.setYRange(1440, 1480);
+        gnuplot_.setXlabel("pressure [Pa]");
+        gnuplot_.setYlabel("density [kg/m^3]");
+        gnuplot_.addDataSetToPlot(x, y, "YourComponent_density.dat", "w l t 'Liquid density of your component'");
+        gnuplot_.plot("YourComponent_density");
+    }
+
     Scalar eps_; //! small epsilon value
     Scalar depthBOR_; //! depth at the bottom of the reservoir
+    Dumux::GnuplotInterface<double> gnuplot_; //! collects data for plotting
 };
 }
 
diff --git a/exercises/exercise-fluidsystem/exercise-fluidsystem_a.input b/exercises/exercise-fluidsystem/exercise-fluidsystem_a.input
index c7c829043c00af4651a245b1b1da7036a6336ad1..8b18ca89524c9fefc63736d232498b67343b6939 100644
--- a/exercises/exercise-fluidsystem/exercise-fluidsystem_a.input
+++ b/exercises/exercise-fluidsystem/exercise-fluidsystem_a.input
@@ -8,3 +8,6 @@ Name = exercise-fluidsystem_a # name will be given to e.g. to the vtk result fil
 [Grid]
 UpperRight = 60 60 # x-/y-coordinates of the upper-right corner of the grid [m]
 Cells = 60 60 # x-/y-resolution of the grid
+
+[Output]
+PlotDensity = false # plot density over pressure for your component