diff --git a/examples/Waste Heat Recovery/heat.py b/examples/Waste Heat Recovery/heat.py
new file mode 100644
index 0000000000000000000000000000000000000000..a02c145ae19f2754af57a3dd3d98cd92cf6fc12a
--- /dev/null
+++ b/examples/Waste Heat Recovery/heat.py	
@@ -0,0 +1,138 @@
+# Purpose: to compute the heat generated by a reaction
+# Mise en stand-by
+
+from thermopy3 import nasa9polynomials as nasa9
+import pandas as pd
+
+temperature = 313.15 #Kelvin
+
+db = nasa9.Database()
+
+#global parameters 
+
+#chemicals 
+H2 = db.set_compound("H2")
+O2 = db.set_compound("O2")
+CH4 = db.set_compound("CH4")
+CO2 = db.set_compound("CO2")
+CO = db.set_compound("carbon monoxide")
+H2OL = db.set_compound("H2O(l)")
+H2O = db.set_compound("H2O")
+
+def compute_reaction_energy(temperature,reactants,products,reactants_coeff,products_coeff):
+    """
+    PRE : 
+        temperature : temperature of the reaction
+        reactants : list of reactants strings matching nasa9 database
+        products : list of products strings matching nasa9 database
+        reactants_coeff : list of coefficients of reactants
+        products_coeff : list of coefficients of products
+    POST : returns the reaction's enthalpy and entropy*T [kJ/mol]
+
+    WARNING : 
+        - components states must be specified, by default set to gas so put H2O(l) for liquid water as H2O stands for water vapor
+        - gibs integrated function of thermopy3 is not working properly, so the function is computed manually
+    """
+    reactants_compounds = [db.set_compound(reactant) for reactant in reactants]
+    products_compounds = [db.set_compound(products) for products in products]
+    # reaction = nasa9.Reaction(temperature, reactants_compounds, products_compounds, reactants_coeff, products_coeff)
+    # print(reaction.gibbs_energy_difference(temperature))
+    deltaH =0
+    deltaS =0
+    for p_idx in range(len(products)):
+        deltaH += products_coeff[p_idx]*products_compounds[p_idx].enthalpy_of_formation
+        deltaS += products_coeff[p_idx]*products_compounds[p_idx].entropy(temperature)
+
+    for r_idx in range(len(reactants)):
+        deltaH -= reactants_coeff[r_idx]*reactants_compounds[r_idx].enthalpy_of_formation
+        deltaS -= reactants_coeff[r_idx]*reactants_compounds[r_idx].entropy(temperature)
+    
+    return  [deltaH/1000 ,deltaS*temperature/1000]
+
+def compute_process_heat(temperature,process_reactants,process_products,process_reactants_coeff,process_products_coeff):
+    """
+    PRE : 
+        temperature : temperature of the reaction
+        process_reactants : list of list containing the reactions'reactants of the global as strings
+        process_products : list of list containing the reactions'products of the global as strings
+        process_reactants_coeff : list of list containing the reactions'reactants' coefficients 
+        process_products_coeff : list of list containing the reactions'products' coefficients
+    POST : returns the heat generated by the reaction [kJ/mol] using compute_reaction_energy()
+
+    WARNING : error while looking for CO!
+    """
+    heat = 0
+    for i in range(len(process_reactants)):
+        heat += compute_reaction_energy(temperature,process_reactants[i],process_products[i],process_reactants_coeff[i],process_products_coeff[i])
+
+def compute_cooling_power_methanation(process_temperature,target_temperature,conversion_efficiency):
+    """
+    PRE : process_temperature : temperature of the methanation process [K]
+          target_temperature : temperature of the methanation's products [K]
+          conversion_efficiency : conversion_efficiency of the methanation process 
+    POST: computes thermal power retrieved from methenation's products cooled down to target temperature [kJ/mol]
+    """
+    reactants = [H2,CO2]
+    products = [CH4,H2OL,H2]
+    total_cooling_heat =0
+    for t in range(target_temperature,process_temperature):
+        total_cooling_heat+=CH4.heat_capacity(t)
+    return total_cooling_heat/1000
+
+def tests():
+
+    
+    print("H2 enthalpy " + str(H2.enthalpy_of_formation) + " Entropy "+ str(H2.entropy(temperature)))
+    print("O2 enthalpy " + str(O2.enthalpy_of_formation) + " Entropy "+ str(O2.entropy(temperature)))
+    print("H2O(l) enthalpy " + str(H2OL.enthalpy_of_formation) + " Entropy "+ str(H2OL.entropy(temperature)))
+
+    print("CH4 enthalpy " + str(CH4.enthalpy_of_formation) + " Entropy "+ str(CH4.entropy(temperature)))
+    print("CO2 enthalpy " + str(CO2.enthalpy_of_formation) + " Entropy "+ str(CO2.entropy(temperature)))    
+    print("CO enthalpy " + str(CO.enthalpy_of_formation) + " Entropy "+ str(CO.entropy(temperature)))
+
+def compute_heat_power(heat_per_mol,mass_output_flux,molar_mass):
+    """
+    PRE : heat_per_mol : heat generated by the reactions inside a technology [kJ/mol]
+          mass_output_flux : mass flux of the technology's output [kt/h]
+          molar_mass : molar mass of the technology's output [g/mol]
+    POST : returns the heat power consummed by the technology [GW]
+    """
+    return heat_per_mol*mass_output_flux*1e9/molar_mass/3600/1e6 #[GW]
+
+def get_data_from_row_col(df,row_name,col_name):
+    """
+    PRE : df : dataframe containing the data
+          row_name : name of the row containing the data
+          col_name : name of the column containing the data
+    POST : returns the data contained in the cell [row_name,col_name] of df
+    """
+    return df.loc[row_name,col_name]
+
+
+def main():
+    
+    #data - Mathias Berger's hub
+    df = pd.read_csv("results/stats/statistics.csv")
+    df.index = df["Unnamed: 0"] #set the index to the first column
+    
+    electrolysis_h2_flux_out = df.loc["ELECTROLYSIS_PLANTS_hydrogen","mean"] #kt/h
+    methanation_ch4_flux_out = df.loc["METHANATION_PLANTS_methane","mean"] #kt/h
+    print("Electrolysis H2 flux out [kt/h]",electrolysis_h2_flux_out)
+    print("Methanation CH4 flux out [kt/h]",methanation_ch4_flux_out)
+
+    electrolysis_energy =  compute_reaction_energy(temperature,["H2O(l)"], ["H2", "O2"],[1], [1, 0.5])
+    methanation_energy = compute_reaction_energy(temperature,["CO2","H2"], ["CH4","H2O"],[1, 4], [1, 2])
+
+    print("Energy consummed by electrolysis [kJ/mol]",electrolysis_energy[0]-electrolysis_energy[1]) #Gibs energy must be considered to set Amperage and voltage
+    print("Heat consummed by methanation [kJ/mol]",methanation_energy[0] ) #Enthalpy is considered
+
+    electrolysis_heat_power = compute_heat_power(electrolysis_energy[0]-electrolysis_energy[1],electrolysis_h2_flux_out,18)
+    methanation_heat_power = compute_heat_power(methanation_energy[0],methanation_ch4_flux_out,16)
+    print("Electrolysis heat power {} [GW]".format(electrolysis_heat_power))
+    print("Methanation heat power {} [GW]".format(methanation_heat_power))
+    # tests()
+    print("Cooling heat",compute_cooling_power_methanation(600,200,0.9))
+if __name__ == "__main__":
+    main()
+
+