With ODEINT, I am trying to configure the methane conversion over the height of a reactor. However, the figure my script produces is far from the figure from the literature. The value directly jumps to a number close to the maximum conversion. The approach is based on the following research.

To confirm I did the coding right, I set the superficial velocity (Jg), gas holdup (GH), and pressure (P) more or less constant in the range of the results from the research.

Can someone identify what leads to the difference between my figure and the figure from the literature??

My current code looks as followed:

def odes(Initial,l): 
    MC = Initial[0]
    P  = Initial[1] #in bar multiply with 101325 to make Pa 
    JG = Initial[2] 
    GH = Initial[3]
    
    #Fixed constants
    R=8.314
    g=9.81
    
    #variable constants
    feed= 2544 #mol/s
    T=1100 #Temperature
    T=T+273.15
    L=14.3 #Length
    D=3.5 #Diameter
    
    #superficial flow
    JG = 0
    
    #gas_holdup
    GH = -(0.241-0.32)/L
    #GH=0.267
    
    #Pressure
    dPdl  = -(57.8-51.5)/L
    
    #Methane conversion constants
    n=1.0809
    k0=1.4676e10 #or 1.4676e4 research states two numbers
    Ea=284948
    k=k0*np.exp(-Ea/(R*T))
    Kc=(101325/(R*T))*np.exp(13.2714-(91204.6/(R*T)))
    
    dMCdl = ((k/feed)*(((1-MC)/(1+MC)*((P*101325)/R*T))**n)*(1-((4*MC**2)/(Kc*(1-MC**2)))*P*101325/(R*T))*(GH*np.pi*D**2)/4) #methane conversion
    
    return[dMCdl,dPdl,JG,GH]
    
Initial=[0,57.8,0.10992158,0.241]
print(odes(Initial=Initial,l=0))
l=np.linspace(0,14.3,144)     
    
Data=odeint(odes,Initial,l) 
fig, ax1 = plt.subplots() 
ax2 = ax1.twinx() 
ax1.plot(l, Data[:,0],label='Methane conversion', color="b") 
ax1.plot(l, Data[:,3],label='Gas holdup',color="g") 
ax2.plot(l,Data[:,1],label='Pressure',color='r')