diff --git a/examples/supernovae_yields/plot_yields.py b/examples/supernovae_yields/plot_yields.py
new file mode 100644
index 0000000000000000000000000000000000000000..2b77499daf97a57df544324789852992ded8934b
--- /dev/null
+++ b/examples/supernovae_yields/plot_yields.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+from time import strftime, gmtime
+plt.rc('text', usetex=True)
+
+
+N = 1000
+m_min = 0.08  # in solar mass
+m_max = 50
+m_wd = 1.4  # white dwarf mass
+solMass_in_cgs = 1.989e33
+
+pyswiftsim.downloadYieldsTable()
+
+
+if __name__ == "__main__":
+    with pyswiftsim.ParameterFile() as filename:
+        parser = pyswiftsim.parseYamlFile(filename)
+
+        # Get masses
+        mass = np.logspace(np.log10(m_min), np.log10(m_max), N,
+                           dtype=np.float32)
+
+        # transform into internal units
+        units = parser["InternalUnitSystem"]
+        unit_mass = units["UnitMass_in_cgs"]
+        unit_vel = units["UnitVelocity_in_cgs"]
+        unit_length = units["UnitLength_in_cgs"]
+        unit_time = unit_length / unit_vel
+
+        mass *= solMass_in_cgs / unit_mass
+
+        m1 = np.ones(N, dtype=np.float32) * 1e-20  # avoid log(0)
+
+        print("Computing supernovae yields...")
+
+        yields_snii = libstellar.supernovaeIIYields(filename, m1, mass)
+        mass_ejected_snii = libstellar.supernovaeIIMassEjected(
+            filename, m1, mass)
+        mass_ejected_np_snii = libstellar.supernovaeIIMassEjectedNonProcessed(
+            filename, m1, mass)
+
+        yields_snia = libstellar.supernovaeIaYields(filename)
+
+        elements = libstellar.getElements(filename)
+
+        # units
+        mass /= solMass_in_cgs / unit_mass
+
+        plt.figure()
+        date = strftime("%d %b %Y", gmtime())
+        plt.title("Date: {}".format(date))
+        for i in range(len(elements)):
+            p = plt.loglog(mass, yields_snii[:, i], label=elements[i])
+            plt.loglog(m_wd, yields_snia[i], "o", color=p[0].get_color())
+
+        plt.loglog(mass, mass_ejected_snii, label="Mass Ejected")
+
+        plt.loglog(mass, mass_ejected_np_snii,
+                   label="Non Processed Mass Ejected")
+
+        plt.xlabel("$Mass [M_\odot]$")
+        plt.ylabel("Mass fraction")
+        plt.legend()
+
+        plt.show()
diff --git a/examples/supernovae_yields/run.sh b/examples/supernovae_yields/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..358d7c530922339afb1903ef7521946d5f6ebeb2
--- /dev/null
+++ b/examples/supernovae_yields/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python plot_yields.py
diff --git a/tests/test_yields.py b/tests/test_yields.py
new file mode 100644
index 0000000000000000000000000000000000000000..f24b4397506f882e04c65fbff984efb32b808f46
--- /dev/null
+++ b/tests/test_yields.py
@@ -0,0 +1,77 @@
+#!/usr/bin/env python3
+# ########################################################################
+# This file is part of PYSWIFT.
+# Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+# ########################################################################
+
+from pyswiftsim import libstellar
+import pyswiftsim
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.rc('text', usetex=True)
+
+
+N = 10
+m_min = 1  # in solar mass
+m_max = 50
+solMass_in_cgs = 1.989e33
+
+pyswiftsim.downloadYieldsTable()
+
+
+if __name__ == "__main__":
+    with pyswiftsim.ParameterFile() as filename:
+        parser = pyswiftsim.parseYamlFile(filename)
+
+        # Get masses
+        mass = np.logspace(np.log10(m_min), np.log10(m_max), N + 1,
+                           dtype=np.float32)
+
+        # transform into internal units
+        units = parser["InternalUnitSystem"]
+        unit_mass = units["UnitMass_in_cgs"]
+        unit_vel = units["UnitVelocity_in_cgs"]
+        unit_length = units["UnitLength_in_cgs"]
+        unit_time = unit_length / unit_vel
+
+        mass *= solMass_in_cgs / unit_mass
+
+        m1 = mass[:-1]
+        m2 = mass[1:]
+
+        print("Computing supernovae yields...")
+
+        yields_snii = libstellar.supernovaeIIYields(filename, m1, m2)
+        mass_ejected_snii = libstellar.supernovaeIIMassEjected(
+            filename, m1, m2)
+        mass_ejected_np_snii = libstellar.supernovaeIIMassEjectedNonProcessed(
+            filename, m1, m2)
+
+        m1 /= solMass_in_cgs / unit_mass
+
+        yields_snia = libstellar.supernovaeIaYields(filename)
+
+        print("Mass used: {}".format(m1))
+
+        print("SNII yields obtained: {}".format(yields_snii))
+
+        print("SNIa yields obtained: {}".format(yields_snia))
+
+        print("SNII mass ejected: {}".format(mass_ejected_snii))
+
+        print("SNII non processed mass ejected: {}".format(
+            mass_ejected_np_snii))