diff --git a/examples/example00.py b/examples/example00.py
new file mode 100644
index 0000000..8f8dfd0
--- /dev/null
+++ b/examples/example00.py
@@ -0,0 +1,72 @@
+#!/usr/bin/env python
+# coding: utf-8
+#
+# Copyright © 2022 - Rennes Physics Institute
+#
+# This file is part of MsSpec-DFM.
+#
+# MsSpec-DFM is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+# MsSpec-DFM 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 General Public License
+# along with MsSpec-DFM. If not, see .
+#
+# Source file : src/python/msspec_dfm/version.py
+# Last modified: Fri, 25 Feb 2022 17:27:32 +0100
+# Committed by : Sylvain Tricot 1645806435 +0100
+
+
+
+# This example shows how to compute the plasmon dispersion for
+# different values of electron densities
+
+# Import section
+from matplotlib import pyplot as plt
+import numpy as np
+
+from msspec_dfm import compute
+
+
+# Electron densities (electrons/cm**3)
+densities = np.array([1e17, 1e18, 1e20, 1e23])
+# Bohr radius (cm)
+a0 = 0.529e-8
+# Temperature
+T = 300
+
+# Loop over density values
+for d in densities:
+ # The density is controlled by the 'RS' parameter, so
+ # compute the corresponding averaged radius (in Bohr units)
+ rs = ( (1/d) / ((4/3) * np.pi) ) ** (1/3) / a0
+ # Some output printing
+ print("Computing {} for RS={}".format(plot_type, rs))
+ # Here is the actual call to calculate
+ data = compute(N_Q=100, RS=rs, T=T)
+ # 'data' is a dict, each key is the name of an output file (without extansion).
+ # The value is the array of output values. The size of this array depends on
+ # the number of points chosen and on the kind of output file.
+
+ # Reading data
+ X, Y = data['plas_disp'].T
+ # Plotting
+ plt.plot(X, Y, label=r'density={:.2e} $e/cm^3$'.format(d))
+
+# Plot setup
+plt.xlabel(r'$q/k_F$')
+plt.ylabel(r'$E/E_F$')
+plt.title('Plasmon Energy [T = {:.1f}K]'.format(T))
+plt.grid(True)
+plt.legend()
+# Pops up the graph window
+plt.show()
+
+
+
+
+