diff --git a/constants_lx_rk1.py b/constants_lx_rk1.py
new file mode 100644
index 0000000000000000000000000000000000000000..8687dd826dc73d3233a7cde645fb598bec743c0d
--- /dev/null
+++ b/constants_lx_rk1.py
@@ -0,0 +1,63 @@
+# -*- coding: utf-8 -*-
+#________________________________Constant values used by Ling Xu_________________________________
+
+#Algorithm constants
+delta_t = 3e-2 #s. It is the time step of the algorithm.
+tf = 500 #s. It is the simulation time.
+
+#Specific constants of the problem
+Tfc = 343 #K. It is the temperature of the fuel cell.
+Pfc = 101325 #Pa. It is the pressure of the fuel cell.
+k = 3 #. It is the ratio of the inlet molar flow of O2 over the reaction molar flow of O2.
+Re = 1e-5 #ohm.m². It is the electron conduction resistance of the circuit.
+
+#The physical characteristics of the PEMFC
+ #Gas channel
+Lch = 1.298 #m. It is the length of the channel.
+Hgc= 1e-3 #m. It is the thickness of the gas channel.
+ #Gas diffusion layer
+Hgdl= 2.1e-4 #m. It is the thickness of the gas diffusion layer.
+epsilon_gdl = 0.6 #It is the anode/cathode GDL porosity.
+ #Catalyst layer
+Hcl = 1e-5 #m. It is the thickness of the anode or cathode catalyst layer.
+epsilon_cl = 0.6 #It is the porosity of the catalyst layer, without units.
+i0_c = 0.1 #A.m-2.It is the exchange current density at the cathode.
+s_lim = 0.3 #It is the limit liquid saturation.
+ #Membrane
+Hmem= 2.5e-5 #m. It is the thickness of the membrane.
+rho_mem = 1980 #kg.m-3. It is the density of the dry membrane.
+M_eq = 1.1e0 #kg.mol-1. It is the equivalent molar mass of ionomer.
+
+#The constants based on the interaction between the structure and the flow
+ #Vapor
+gamma_v = 1.3 #s-1. It is the sorption/desorption rate of vapor at the ionomer/CL interface.
+Dvc = 2.236e-5 #m².s-1. It is the vapor diffusivity in cathode.
+Dva = 5.457e-5 #m².s-1. It is the vapor diffusivity in anode.
+h_mv = 0.08 #m.s-1. It is the convective mass transfer coefficient of vapor.
+ #Liquid water
+gamma_l = 100 #s-1. It is the sorption/desorption rate of liquid at the ionomer/CL interface.
+Kc = 6.875e-13 #m². It is the anode/cathode GDL permeability for liquid water.
+theta_c = 110 #°. It is the contact angle of GDL for liquid water.
+ #Dioxygen
+D_O2 = 2.9e-5 #m².s-1. It is the O2 diffusivity in cathode.
+h_O2 = 0.078 #m.s-1. It is the convective mass transfer coefficient of O2.
+
+#The physical constants
+F = 96485 #C.mol-1. It is the Faraday constant.
+R = 8.314 #J. mol-1.K-1. It is the universal gas constant.
+M_H2O = 0.018 #kg.mol-1. It is the water molecular mass.
+rho_H2O = 1000 #kg.m-3. It is the water density.
+nu_l = 3.7e-7 #m².s-1. It is the liquid water kinematic viscosity.
+mu_l = 3.56e-4 #Pa.s. It is the liquid water dynamic viscosity.
+sigma = 0.0625 #N.m-1. It is the water surface tension.
+Csat = 10.9 #mol.m-3. It is the saturated vapor concentration at 70°C for a perfect gas.
+Psat = 31200 #Pa. It is the saturated partial pressure of vapor at 70°C.
+C_O2ref = 40.88 #mol.m-3. It is the reference concentration of oxygen.
+y_O2 = 0.2095 #. It is the molar fraction of O2 in dry air.
+
+#Mathematical factors
+Kshape = 5 #It is a shape mathematical factor, without physical sense.
+
+#Not understood yet
+mu_cg = 1.881e-5 #Pa.s. It is the cathode gas dynamic viscosity.
+alpha_c = 0.5 #It is the transfer coefficient of the cathode.
\ No newline at end of file
diff --git a/current_density.py b/current_density.py
new file mode 100644
index 0000000000000000000000000000000000000000..4c8af8779258aae3ae2f7f826991218c9330d051
--- /dev/null
+++ b/current_density.py
@@ -0,0 +1,21 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+#_______________________________________Current density function______________________________________
+
+#The current density raises in ~5s from 0.0 to 0.4 A.m-2 at tf/2.
+
+def current_density(n,delta_t):
+
+ i_fc = np.zeros(n)
+ c = round(n/2) #The current density value changes at half-time.
+ tlim=c*delta_t #time at which the current density raises.
+ tau_c=1 #In 3*tau_c, 95% of the final value is reached.
+
+ for j in range(c):
+ i_fc[j]=0.0 #A.m-2. Start at 0.0 A.m-2
+
+ for j in range(c,n):
+ t=j*delta_t #current time.
+ i_fc[j]=0.0+0.4*(1-np.exp( -(t-tlim)/tau_c )) #A.m-2.
+
+ return i_fc
\ No newline at end of file
diff --git a/main_rk1.py b/main_rk1.py
new file mode 100644
index 0000000000000000000000000000000000000000..47e050c9cc30ec83e970a6c2f54d36ab4f4d6be8
--- /dev/null
+++ b/main_rk1.py
@@ -0,0 +1,443 @@
+# -*- coding: utf-8 -*-
+from IPython import get_ipython
+import sys
+import time
+import numpy as np
+import numpy.polynomial.polynomial as nppol
+import matplotlib.pyplot as plt
+
+from current_density import current_density
+
+#Clean the console and reset all the variables
+get_ipython().magic('clear')
+get_ipython().magic('reset -f')
+
+#Start time
+start_time = time.time()
+
+#Select the value of the constants needed for the simulation.
+#"lx" : data used by Ling Xu in his articles + other data to complete them.
+version = "lx"
+
+#_______________________________Initialisation and constants value______________________________
+
+if version == "lx":
+ from constants_lx_rk1 import Re,delta_t,tf,Tfc,Pfc,k,Hgc,Hgdl,epsilon_gdl,Hcl,epsilon_cl,i0_c,s_lim,\
+ Hmem,rho_mem,M_eq,gamma_v,Dvc,Dva,h_mv,Kc,theta_c,D_O2,h_O2,F,R,M_H2O,rho_H2O,nu_l,mu_l,\
+ sigma,Csat,Psat,C_O2ref,y_O2,mu_cg,alpha_c
+ from transitory_functions_lx import lambda_eq,a_w,Cw,gamma,D,Q
+
+#Initialisation. Is has to start in the vapor regime ! (for now ?)
+Cini = 10.6
+C_O2ini = 1.55
+s5 = 0
+
+Cagc,C1,C5,Ccgc,C_O2cgc = Cini,Cini,Cini,Cini,C_O2ini
+lambda2,lambda3,lambda4 = lambda_eq(a_w(Cini)),lambda_eq(a_w(Cini)),lambda_eq(a_w(Cini))
+
+#Total number of step
+n=round(tf/delta_t)
+
+#Target variables
+i_fc_t=current_density(n,delta_t) #table with all the values of i_fc
+Ja_t,Jc_t,J_O2_t,Jmem_acl_t,Jmem_ccl_t = np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n)
+Cagc_t,C0_t,C1_t,lambda2_t,lambda3_t = np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n)
+lambda4_t,s5_t,C5_t,C6_t,Ccgc_t = np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n)
+C_O2cgc_t,C_O2inter_t,C_O2ccl_t,Vcell_t = np.zeros(n),np.zeros(n),np.zeros(n),np.zeros(n)
+regime_t,a_w_c_t,lambda_eq_t = [],np.zeros(n),np.zeros(n)
+
+#_______________________________________Loop implementation_______________________________________
+for j in range (0,n):
+
+ #____________________________________Current density value____________________________________
+ i_fc=i_fc_t[j]
+
+ #_______________________Concentration flows functions of i_fc_______________________
+
+ if j==0 or i_fc != i_fc_t[j-1]: #otherwise it is useless to recalculate them
+ #Intermediates
+ Ccgc_eq = Csat - Psat/(Pfc-Psat)*i_fc/(2*F)*( Hgdl/(Dvc*epsilon_gdl**1.5) + 1/h_mv )
+ Phi_eq = Ccgc_eq*R*Tfc/Psat
+
+ #Concentration flows at the GC
+ Wagc_in = Psat/(Pfc-Psat)*i_fc/(2*F*Hgc)
+ W_O2_in = k*i_fc/(4*F*Hgc)
+ W_O2_out = (k-1)*i_fc/(4*F*Hgc)
+ Wcgc_in = Phi_eq*Psat/(Pfc-Phi_eq*Psat)*k*i_fc/(4*F*y_O2*Hgc)
+ Wcgc_out = i_fc/(2*F*Hgc)*(1+ Psat/(Pfc-Psat) + Phi_eq*Psat/(Pfc-Phi_eq*Psat)*k/(2*y_O2))
+
+ #Oxygen transport flow in the CGDL:
+ J_O2 = -i_fc/(4*F)
+
+ #__________________________________Calculation of the flows__________________________________
+
+ #Water transport flow in the membrane
+ Jmem_acl = 2.5/22*i_fc/F*lambda2 - 2*rho_mem/M_eq*D(lambda2)*(lambda3-lambda2)/Hmem
+ Jmem_ccl = 2.5/22*i_fc/F*lambda4 - 2*rho_mem/M_eq*D(lambda4)*(lambda4-lambda3)/Hmem
+
+ #Intermediate values
+ Cw_c = Cw(s5,C5)
+ a_w_a, a_w_c= a_w(C1), a_w(Cw_c)
+ gamma_c = gamma(a_w_c)
+
+ #Water transport flow in the two GDLs:
+ Ja = gamma_v*epsilon_cl*Hcl*rho_mem/M_eq*(lambda_eq(a_w_a)-lambda2)
+
+ if lambda4 >= 22:
+ Jc = max( Jmem_ccl + i_fc/(2*F) ,\
+ gamma_c*epsilon_cl*Hcl*rho_mem/M_eq*(lambda4-lambda_eq(a_w_c)) )
+ else:
+ Jc = gamma_c*epsilon_cl*Hcl*rho_mem/M_eq*(lambda4-lambda_eq(a_w_c))
+
+ #________________________Determination of the water transport regime_______________________
+
+ if C5 < Csat :
+ regime = "V"
+ elif C5 >= Csat and Ccgc < Csat :
+ regime = "M"
+ elif C5 >= Csat and Ccgc >= Csat :
+ regime = "L"
+ else:
+ print(C5,Ccgc,Csat)
+ sys.exit("The concentrations are not relevant")
+
+ #________________________Calculation of the instantaneous quantities________________________
+
+ #At the anode :
+ C0 = Cagc - Ja/h_mv
+ if C0>Csat:
+ C0=Csat
+ if C0<0:
+ print (C0,j)
+ sys.exit("C0 concentration cannot be negative.")
+
+ C1 = C0 - Ja*Hgdl/(Dva*epsilon_gdl**1.5)
+ if C1>Csat:
+ C1=Csat
+ if C1<0:
+ print (C1,j)
+ sys.exit("C1 concentration cannot be negative.")
+
+ #At the cathode :
+ if regime == "V":
+ s5,s6 = 0,0
+
+ C6 = Ccgc + Jc/h_mv
+ if C6>Csat:
+ C6=Csat
+ if C6<0:
+ print (C6,j)
+ sys.exit("C6 concentration cannot be negative.")
+
+ C5 = C6 + Jc*Hgdl/(Dvc*epsilon_gdl**1.5)
+ if C5>Csat:
+ C5=Csat
+ if C5<0:
+ print (C5,j)
+ sys.exit("5 concentration cannot be negative.")
+
+ C_O2inter = C_O2cgc + J_O2/h_O2
+ if C_O2inter<0:
+ print (C_O2inter,j)
+ sys.exit("C_O2inter concentration cannot be negative.")
+
+ C_O2ccl = C_O2inter + J_O2*Hgdl/(D_O2*epsilon_gdl**1.5)
+ if C_O2ccl<0:
+ print (C_O2ccl,j)
+ sys.exit("C_O2ccl concentration cannot be negative.")
+
+
+ elif regime == "M":
+ C5,s6 = Csat,0
+
+ C6 = Ccgc + Jc/h_mv
+ if C6>Csat:
+ C6=Csat
+ if C6<0:
+ print (C6,j)
+ sys.exit("C6 concentration cannot be negative.")
+
+ xs = Hgdl + Dvc*epsilon_gdl**1.5*(C6-Csat)/Jc
+
+ a0 = M_H2O*Jc*xs / (sigma*Kc/nu_l*np.cos(theta_c)*(epsilon_gdl/Kc)**0.5)
+ solv = nppol.polyroots([a0, 0, 0, 0, 0.35425, -0.848, 0.6315])
+ c=0; #number of time it enters in the following if :
+ for i in solv:
+ if abs(i.imag) < 1e-10 and i.real>=0.0 and i.real<=1.0001:
+ if s5 >= 1.0:
+ s5=0.999999
+ else :
+ s5 = i.real
+ c += 1
+ if c != 1 :
+ print (solv,j)
+ sys.exit("The polynomial solver gives more or less than a unique useful value.")
+
+ C_O2inter = C_O2cgc + J_O2/h_O2
+ if C_O2inter<0:
+ print (C_O2inter,j)
+ sys.exit("C_O2inter concentration cannot be negative.")
+
+ C_O2_xs = C_O2inter + J_O2*(Hgdl-xs)/(D_O2*epsilon_gdl**1.5)
+ if C_O2_xs<0:
+ print (C_O2_xs,j)
+ sys.exit("C_O2_xs concentration cannot be negative.")
+
+ C_O2ccl = C_O2_xs - (J_O2*sigma*Kc/nu_l*np.cos(theta_c)*(epsilon_gdl/Kc)**0.5) / \
+ (M_H2O*Jc*D_O2*epsilon_gdl**1.5) * (Q(s5)-Q(0))
+ if C_O2ccl<0:
+ print (C_O2ccl,j)
+ sys.exit("C_O2ccl concentration cannot be negative.")
+
+ elif regime == "L":
+ C5,C6 = Csat,Csat
+
+ s6 = 1/(1+ (rho_H2O*Wcgc_out*Hgc*mu_cg / (M_H2O*abs(Jc)*Ccgc_eq*mu_l))**(1/3) )
+
+ a0 = M_H2O*Jc*Hgdl / (sigma*Kc/nu_l*np.cos(theta_c)*(epsilon_gdl/Kc)**0.5) \
+ -0.35425*s6**4 + 0.848*s6**5 -0.6315*s6**6
+ solv = nppol.polyroots([a0, 0, 0, 0, 0.35425, -0.848, 0.6315])
+ c=0; #number of time it enters in the following if :
+ for i in solv:
+ if abs(i.imag) < 1e-10 and i.real>=0.0 and i.real<=1.0001:
+ if s5 >= 1.0:
+ s5=0.999999
+ else :
+ s5 = i.real
+ c += 1
+ if c != 1 :
+ print (solv,j)
+ sys.exit("The polynomial solver gives more or less than a unique useful value.")
+
+ C_O2inter = C_O2cgc + J_O2/h_O2
+ if C_O2inter<0:
+ print (C_O2inter,j)
+ sys.exit("C_O2inter concentration cannot be negative.")
+
+ C_O2ccl = C_O2inter - (J_O2*sigma*Kc/nu_l*np.cos(theta_c)*(epsilon_gdl/Kc)**0.5) / \
+ (M_H2O*Jc*D_O2*epsilon_gdl**1.5) * (Q(s5)-Q(s6))
+ if C_O2ccl<0:
+ print (C_O2ccl,j)
+ sys.exit("C_O2ccl concentration cannot be negative.")
+
+ #___________________________Control of the water transport regime___________________________
+
+ #I am going to look at this later.
+ # if regime == "V":
+ # if C5 > Csat :
+ # regime == "M"
+ # elif regime == "M":
+ # if C5 < Csat:
+ # regime == "V"
+ # if C5 >= Csat and Ccgc >= Csat :
+ # regime == "L"
+ # elif regime == "L":
+ # if Ccgc < Csat:
+ # regime == "M"
+
+ #_____________________________Calculation of the slow quantities_____________________________
+
+ Cagc += delta_t*(Wagc_in - Ja/Hgc)
+ if Cagc>Csat:
+ Cagc=Csat
+ if Cagc<0:
+ print (Cagc,j)
+ sys.exit("Cagc concentration cannot be negative.")
+
+ Ccgc += delta_t*(Wcgc_in + Jc/Hgc - Wcgc_out)
+ if Ccgc>Csat:
+ Ccgc=Csat
+ if Ccgc<0:
+ print (Ccgc,j)
+ sys.exit("Ccgc concentration cannot be negative.")
+
+ C_O2cgc += delta_t*(W_O2_in + J_O2/Hgc - W_O2_out)
+ if C_O2cgc<0:
+ print (C_O2cgc,j)
+ sys.exit("C_O2cgc concentration cannot be negative.")
+
+ lambda2 += delta_t*M_eq/(Hcl*rho_mem)*(Ja-Jmem_acl)
+ if lambda2<0:
+ print (lambda2,j)
+ sys.exit("lambda2 concentration cannot be negative.")
+
+ lambda3 += delta_t*M_eq/(Hmem*rho_mem)*(Jmem_acl-Jmem_ccl)
+ if lambda3<0:
+ print (lambda3,j)
+ sys.exit("lambda3 concentration cannot be negative.")
+
+ lambda4 += delta_t*M_eq/(Hcl*rho_mem)*(Jmem_ccl-Jc+i_fc/(2*F))
+ if lambda4<0:
+ print (lambda4,j)
+ sys.exit("lambda4 concentration cannot be negative.")
+
+ #______________________________Calculation of the cell voltage______________________________
+
+ Ueq = 1.229 - 8.5e-4*(Tfc-298.15)
+
+ if i_fc == 0.0:
+ eta_c = 0.0
+ Rccl = 0.0
+ Rmem = 0.0
+ else:
+ eta_c = -R*Tfc/(alpha_c*F)*np.log(s_lim/(s_lim-s5) * C_O2ref/C_O2ccl * i_fc/i0_c )
+ Rccl = Hcl/( epsilon_cl**1.5 * (0.5139*lambda4-0.326) * np.exp(1268*(1/303.15 - 1/Tfc)) )
+ Rmem = Hmem/( 2*np.exp(1268*(1/303.15 - 1/Tfc)) ) *\
+ ( np.log(0.5139*lambda4-0.326)-np.log(0.5139*lambda3-0.326) ) / (0.5139*(lambda4-lambda3))
+
+ Rp = Rmem + Rccl/3
+
+ Vcell = Ueq + eta_c - i_fc*(Rp+Re)
+
+ #_______________________________Record of the target variables_______________________________
+
+ Ja_t[j],Jc_t[j],J_O2_t[j],Jmem_acl_t[j],Jmem_ccl_t[j] = Ja,Jc,J_O2,Jmem_acl,Jmem_ccl
+ Cagc_t[j],C0_t[j],C1_t[j],lambda2_t[j],lambda3_t[j] = Cagc,C0,C1,lambda2,lambda3
+ lambda4_t[j],s5_t[j],C5_t[j],C6_t[j],Ccgc_t[j] = lambda4,s5,C5,C6,Ccgc
+ C_O2cgc_t[j],C_O2inter_t[j],C_O2ccl_t[j],Vcell_t[j] = C_O2cgc,C_O2inter,C_O2ccl,Vcell
+ a_w_c_t[j],lambda_eq_t[j] = a_w_c,lambda_eq(a_w_c)
+ regime_t += [regime]
+
+#_____________________________________________Display_____________________________________________
+
+t = np.linspace(0,tf,n)
+
+#Zoom on the timeline [0.49*tf,0.52*tf]s
+n1_z1=round(0.49*tf/delta_t) #number of step to reach 245s.
+n2_z1=round(0.52*tf/delta_t) #number of step to reach 260s.
+delta_n = n2_z1-n1_z1 #number of step from 245s to reach 260s.
+t_z1 = np.linspace(0.49*tf,0.52*tf,delta_n)
+
+#Zoom on the timeline [0.4*tf,0.7*tf]s
+n1_z2=round(0.4*tf/delta_t) #number of step to reach 245s.
+n2_z2=round(0.7*tf/delta_t) #number of step to reach 260s.
+delta_n = n2_z2-n1_z2 #number of step from 245s to reach 260s.
+t_z2 = np.linspace(0.4*tf,0.7*tf,delta_n)
+
+#Zoom on the timeline [0.0*tf,0.01*tf]s
+n1_z3=round(0.0*tf/delta_t) #number of step to reach 245s.
+n2_z3=round(0.01*tf/delta_t) #number of step to reach 260s.
+delta_n = n2_z3-n1_z3 #number of step from 245s to reach 260s.
+t_z3 = np.linspace(0.0*tf,0.01*tf,delta_n)
+
+# #Test
+# n1_z4=n//2
+# n2_z4=25038
+# delta_n = n2_z4-n1_z4 #number of step from 245s to reach 260s.
+# t_z4 = np.linspace(n1_z4*delta_t,n2_z4*delta_t,delta_n)
+
+# plt.figure(80)
+# plt.plot(t_z4,a_w_c_t[n1_z4:n2_z4],label='a_w_c')
+# plt.legend()
+
+# plt.figure(81)
+# plt.plot(t_z4,lambda_eq_t[n1_z4:n2_z4],label='lambda_eq')
+# plt.legend()
+
+# #Figure 1 : zoom1 on i_fc
+# plt.figure(1)
+# i_fc_z=i_fc_t[n1_z1:n2_z1] #It is a zoom on i_fc.
+# plt.plot(t_z1,i_fc_z,label='i_fc')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Current density (A/m²)')
+# plt.title('The current density\nbehaviour over time')
+# plt.legend()
+
+#Figure 2 : the different flows J
+plt.figure(2)
+i_fc_flow = i_fc_t/(2*F)
+plt.plot(t_z1,Ja_t[n1_z1:n2_z1],label='Ja')
+plt.plot(t_z1,Jc_t[n1_z1:n2_z1],label='Jc')
+plt.plot(t_z1,i_fc_flow[n1_z1:n2_z1],label='i_fc/2F')
+plt.plot(t_z1,Jmem_acl_t[n1_z1:n2_z1],label='Jmem_acl')
+plt.plot(t_z1,Jmem_ccl_t[n1_z1:n2_z1],label='Jmem_ccl')
+# plt.plot(t,J_O2_t,label='J_O2')
+plt.xlabel('Time (s)')
+plt.ylabel('Flow (mol/m²/s)')
+plt.title('The flow\nbehaviour over time')
+plt.legend()
+
+# #Figure 3 : the different vapor concentration C
+# plt.figure(3)
+# plt.plot(t,Cagc_t,label='Cagc')
+# plt.plot(t,C0_t,label='C0')
+# plt.plot(t,C1_t,label='C1')
+# plt.plot(t,C5_t,label='C5')
+# plt.plot(t,C6_t,label='C6')
+# plt.plot(t,Ccgc_t,label='Ccgc')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Vapor concentration (mol/m^3)')
+# plt.title('The vapor concentration\nbehaviour over time')
+# plt.legend()
+
+# #Figure 4 : the water content
+# plt.figure(4)
+# plt.plot(t,lambda2_t,label='lambda2')
+# plt.plot(t,lambda3_t,label='lambda3')
+# plt.plot(t,lambda4_t,label='lambda4')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Water content')
+# plt.title('The three water content\nbehaviour over time')
+# plt.legend()
+
+# # Figure 5 : the liquid water saturation s5
+# plt.figure(5)
+# plt.plot(t,s5_t,label='s5')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Liquid water saturation at CCL')
+# plt.title('The CCL liquid water saturation\nbehaviour over time')
+# plt.legend()
+
+# #Figure 6 : the O2 concentrations C_O2
+# plt.figure(6)
+# plt.plot(t,C_O2cgc_t,label='C_O2cgc')
+# plt.plot(t,C_O2inter_t,label='C_O2inter')
+# plt.plot(t,C_O2ccl_t,label='C_O2ccl')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Oxygen concentration (mol/m3)')
+# plt.title('The oxygen concentrations\nbehaviour over time')
+# plt.legend()
+
+# # Figure 7 : zoom1 on Vcell
+# plt.figure(7)
+# Vcell_z=Vcell_t[n1_z1:n2_z1] #It is a zoom on Vcell.
+# plt.plot(t_z1,Vcell_z,label='Vcell')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Cell voltage (V)')
+# plt.title('The cell voltage\nbehaviour over time')
+# plt.legend()
+
+# #Figure 8 : the different vapor concentration C zoom1
+# plt.figure(8)
+# plt.plot(t_z1,Cagc_t[n1_z1:n2_z1],label='Cagc')
+# plt.plot(t_z1,C0_t[n1_z1:n2_z1],label='C0')
+# plt.plot(t_z1,C1_t[n1_z1:n2_z1],label='C1')
+# plt.plot(t_z1,C5_t[n1_z1:n2_z1],label='C5')
+# plt.plot(t_z1,C6_t[n1_z1:n2_z1],label='C6')
+# plt.plot(t_z1,Ccgc_t[n1_z1:n2_z1],label='Ccgc')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Vapor concentration (mol/m^3)')
+# plt.title('The vapor concentration\nbehaviour over time')
+# plt.legend()
+
+# #Figure 9 : the water content zoom 1
+# plt.figure(9)
+# plt.plot(t_z1,lambda2_t[n1_z1:n2_z1],label='lambda2')
+# plt.plot(t_z1,lambda3_t[n1_z1:n2_z1],label='lambda3')
+# plt.plot(t_z1,lambda4_t[n1_z1:n2_z1],label='lambda4')
+# # plt.plot(t_z1,lambda_eq(a_w_c_t)[n1_z1:n2_z1],label='lambda_eq(a_w_c)')
+# plt.xlabel('Time (s)')
+# plt.ylabel('Water content')
+# plt.title('The three water content\nbehaviour over time')
+# plt.legend()
+
+# #Figure 10 : test
+# plt.figure(10)
+# test=C5_t
+# plt.plot(t_z1,test[n1_z1:n2_z1],label='C5')
+# plt.xlabel('Time (s)')
+# plt.legend()
+
+#Algorithm time
+algo_time = time.time() - start_time
+print('Time of the algorithm in second :',algo_time)
\ No newline at end of file
diff --git a/transitory_functions_lx.py b/transitory_functions_lx.py
new file mode 100644
index 0000000000000000000000000000000000000000..17e00c5900d049699fc5eaa020fe6a17c44f0382
--- /dev/null
+++ b/transitory_functions_lx.py
@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+#_________________________________Initialisation and constants value________________________________
+from constants_lx_rk1 import Tfc,gamma_v,gamma_l,M_H2O,rho_H2O,Csat,Kshape
+
+#_______________________________________Transitory functions_______________________________________
+
+def lambda_eq(a_w):
+ return 0.5*( 0.043 +17.81*a_w -39.85*a_w**2 +36.0*a_w**3 )*( 1 -np.tanh(100*(a_w -1)) ) \
+ + 0.5*( 14 +8*( 1-np.exp(-Kshape*(a_w -1)) ) )*( 1+np.tanh(100*(a_w -1)) )
+
+def a_w(Cw):
+ return Cw/Csat
+
+def Cw(s,C):
+ return C+s*rho_H2O/M_H2O
+
+def gamma(a_w):
+ ome = omega(a_w)
+ if a_w <=1 :
+ return gamma_v
+ else:
+ return (1-ome)*gamma_v + ome*gamma_l
+
+def omega(a_w):
+ lambdaEq = lambda_eq(a_w)
+ return (lambdaEq-14)/(22-14)*0.5 + (a_w -1)/(3 -1)*0.5
+
+def D(lambdaa): #the diffusion coefficient of water in the membrane
+ if lambdaa <= 3:
+ return 3.1e-7*lambdaa*(np.exp(0.28*lambdaa)-1)*np.exp(-2346/Tfc)
+ else:
+ return 4.17e-8*lambdaa*(161*np.exp(-lambdaa)+1)*np.exp(-2346/Tfc)
+
+def Q(s):
+ return -(14735*s**5-300*s**4+9439*s**3+18878*s**2+75512*s-151024) / (17500*np.sqrt(1-s))
\ No newline at end of file