Need help with troubleshooting and how can I fix these errors in my code?

Good day, I am trying to emulate the works of Catalan and Thanh Le in their modelling of a methane pyrolysis reactor. I have attached their relevant papers should anyone be interested to read. The m files which contains the code to open in Matlab is also included. The code I came up with is also provided below and is March17Code.m in the attachment. I am trying to use ode45 here but it is not working.

% Define paramters

T = 1150; % Temperature [K]

D = 0.03; % Diameter [m]

P_initial = 101325 * 3; % Initial pressure (Pa)

L_span = [0 2]; % Discretize L from 0 to 1

literFlowRate = 0.5; % flowrate in L/min

epsilon_CH4 = 0.8; % mole fraction of CH4 in reactor feed

epsilon_H2 = 0; % mole fraction of CH4 in reactor feed

epsilon_Ni = 0.27; % mole fraction of Ni in molten bath

V_gj = 0.1; % Drift velocity if j_g < 0.5

% Define constants

R = 8.314; % Universal gas constant (J/(mol·K))

g = 9.81; % Acceleration due to gravity (m/s^2)

A = -1.39; % Empirical constant

No = 6.02214 * 10^23; % Avogadro number

M_Ni = 58.63/1000; % Molar mass of Ni (kg/mol)

M_Bi = 208.98/1000; % Molar mass of Bi (kg/mol)

M_CH4 = 16.04/1000; % Molar mass of CH4 (kg/mol)

M_H2 = 2.02/1000; % Molar mass of H2 (kg/mol)

M_I = 39.95/1000; % Molar mass of Ar assuming Ar is I (kg/mol)

T_Melting = 1013; % Melting temperature of metal (K)

tolerance = 1e-6; % Convergence tolerance

max_iterations = 1000; % Maximum number of iterations

k_nf_0 = 14676000000; % Non-catalytic rate constant

E_nf_a = 284948; % Activation energy in J/mol

k_cf_0 = 78800; % Catalytic rate constant

E_cf_a = 208000; % Activation energy in J/mol

K = exp(13.2714 - (91204.6 / (R * T))); % Equilibrium constant

P0 = 101325; % Pressure

K_c = K * (P0 / (R * T)); % Adjust equilibrium constant

n = 1.0809;

% Initialize parameters

P = P_initial; % Initialize P

X_CH4_initial = 0; % Initialize X_CH4

alphaInitialGuess = 0.5; % Initial guess for alpha_g

% Calculated parameters

n_total = (literFlowRate / 22.4) / 60; % convert L/min to mole/s

n_CH4_b = epsilon_CH4 * n_total;

n_H2_b = epsilon_H2 * n_total;

n_I_b = (1 - epsilon_CH4 - epsilon_H2) * n_total;

totalMoleFlowRateAtBottom = n_CH4_b + n_H2_b + n_I_b;

fraction_H2_b = n_H2_b / n_CH4_b;

fraction_I_b = n_I_b / n_CH4_b;

rho_Ni = 9091 - 0.4932 * T;

rho_Bi = 10698 - 1.2064 * T;

k_cf = k_cf_0 * exp(-E_cf_a/(R * T));

k_nf = k_nf_0 * exp(-E_nf_a/(R * T));

% Density of liquid metal calculation

numerator_rho_L = epsilon_Ni * M_Ni + (1 - epsilon_Ni) * M_Bi;

denominator_rho_L = ((epsilon_Ni * M_Ni) / rho_Ni) + (((1 - epsilon_Ni) * M_Bi) / rho_Bi);

rho_L = numerator_rho_L / denominator_rho_L; % Liquid density (kg/m^3)

% Surface tension of liquid metal calculation

sigma_Bi = 0.4208 - 8.1e-5 * T;

V_Bi = M_Bi / rho_Bi;

A_Bi = 1.091 * (No^(1/3)) * (V_Bi^(2/3));

sigma_L = (sigma_Bi + ((R * T) / A_Bi) * log((1.21 * 0.57)/(1.33*0.43))); % Surface tension (N/m)

% Viscosity of liquid metal calculation

parameterB = 2.65 * (T_Melting^(1.27));

M_L = (epsilon_Ni * M_Ni) + ((1 - epsilon_Ni) * M_Bi);

numerator_parameterA = 1.7e-7 * (rho_L^(2/3)) * (T_Melting^(1/2)) * M_L^(-1/6);

denominator_parameterA = exp( parameterB / (R * T_Melting));

parameterA = numerator_parameterA / denominator_parameterA;

mew_L = parameterA * exp ( parameterB / (R * T)); % Liquid dynamic viscosity (Pa·s)

nu_L = mew_L / rho_L; % kinematic viscosity (m^2/s)

[L, X_CH4_profile] = ode45(@(L, X_CH4) ...

odefunction(L, X_CH4, totalMoleFlowRateAtBottom, k_cf, k_nf, K_c, ...

P_initial, epsilon_CH4, fraction_H2_b, fraction_I_b, R, T, D, rho_L, ...

sigma_L, g, nu_L), L_span, X_CH4_initial);

% Loop over each point in L

function dX_CH4_dL = odefunction(L, X_CH4, totalMoleFlowRateAtBottom, k_cf, k_nf, K_c, ...

epsilon_CH4, fraction_H2_b, fraction_I_b, R, T, D, rho_L, sigma_L, g, nu_L)

P = P_initial;

% Superficial gas velocity

j_g = ((4 * totalMoleFlowRateAtBottom * (1 + epsilon_CH4 * X_CH4)) / (pi * (D^2)))* ((R*T) / P);

% Gas density calculation

rho_g_Numerator = ((fraction_H2_b + 2 * X_CH4) * M_H2) + ((1-X_CH4) * M_CH4) + (fraction_I_b * M_I);

rho_g_Denominator = 1 + fraction_H2_b + X_CH4 + fraction_I_b;

rho_g = (rho_g_Numerator / rho_g_Denominator) * (P / (R * T));

% Gas holdup calculation

j_g_plusSubNumerator = sigma_L * g * abs(rho_L - rho_g);

j_g_plusSubDenominator = rho_L^2;

j_g_plusDenominator = j_g_plusSubNumerator / j_g_plusSubDenominator;

j_g_plus = j_g / ((j_g_plusDenominator^(0.25)));

if j_g_plus < 0.5

V_gj_plus2_SubNumerator = sigma_L * g * abs(rho_L - rho_g);

V_gj_plus2_SubDenominator = rho_L^2;

V_gj_plus2_Denominator = V_gj_plus2_SubNumerator / V_gj_plus2_SubDenominator;

V_gj_plus2 = V_gj / ((V_gj_plus2_Denominator^(0.25)));

elseif j_g_plus > 0.5

N_mewDenominator = rho_L * sigma_L * ((sigma_L/(g * abs(rho_L - rho_g)))^(0.5));

N_mew = nu_L / N_mewDenominator;

D_H_starDenominator = (sigma_L / (g * abs(rho_L - rho_g)));

D_H_star = D / D_H_starDenominator;

if N_mew <= 2.2e-3 && D_H_star <= 30

V_gj_plus2 = 0.0019 * (D_H_star*0.809) * ((rho_g/rho_L)^(-0.157)) * (N_mew^(-0.562));

elseif N_mew <= 2.2e-3 && D_H_star >= 30

V_gj_plus2 = 0.03 * ((rho_g/rho_L)^(-0.157)) * (N_mew^(-0.562));

elseif N_mew >= 2.2e-3 && D_H_star >= 30

V_gj_plus2 = 0.92 * ((rho_g/rho_L)^(-0.157));

end

end

C_0 = 1.2 - (0.2 * ((rho_g / rho_L)^(0.5)));

alpha = alphaInitialGuess; % Initial guess for alpha_g

for iteration = 1:max_iterations

V_gj_plus1 = (2^(0.5)) * ((1- alpha)^1.75);

V_gj_plus = (V_gj_plus1 * exp(A * j_g_plus)) + (V_gj_plus2 * (1 - exp(A * j_g_plus)));

alpha_new = j_g_plus / ((C_0 * j_g_plus) + V_gj_plus);

% Check for convergence

if abs(alpha_new - alpha) < tolerance

break;

end

% Update alpha_g for next iteration

alpha = alpha_new;

% Stop if maximum iterations are reached

if iteration == max_iterations

warning('Maximum iterations reached at L = . Alpha_g may not have converged.');

end

end

% Specific interfacial area calculation

N_Bo = (g * (D^2) * rho_L) / sigma_L;

N_Ga = (g * (D^3)) / (nu_L^2);

N_Fr = j_g / ((g * D)^(0.5));

d_vs = D * 26 * (N_Bo^(-0.5)) * (N_Ga^(-0.12)) * (N_Fr^(-0.12));

a_g = 6 / d_vs;

% Update pressure at current L

dP_dL = -((rho_L * (1 - alpha)) + (rho_g * alpha)) * g;

P = P + dP_dL * (L(2) - L(1)); % Incremental pressure update

% X_CH4 as a function of L

term0 = (alpha * pi * D^2) / 4;

term1_Numerator = a_g * k_cf * (1 - X_CH4);

term1_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4));

term1 = (term1_Numerator / term1_Denominator) * (P / (R * T));

subTerm2_Numerator = 1 - X_CH4;

subTerm2_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4));

subTerm2 = (subTerm2_Numerator / subTerm2_Denominator) * (P / (R * T));

term2 = k_nf * (subTerm2^n);

subTerm3_Numerator = n_CH4_b * ((fraction_H2_b + 2 * X_CH4)^2);

subTerm3_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4)) * (1 - X_CH4) * K_c;

subTerm3 = (subTerm3_Numerator / subTerm3_Denominator) * (P / (R * T));

term3 = 1 - subTerm3;

dX_CH4_dL = term0*(term1 + term2)*(term3);

end

% Plot CH4 Conversion along Reactor Length

figure;

plot(L, X_CH4_profile, 'b-', 'LineWidth', 2);

xlabel('Reactor Length (L) [m]');

ylabel('CH_4 Conversion (X_{CH4})');

title('CH_4 Conversion Profile along Reactor Length');

grid on;

legend('X_{CH4} vs L');

When I try to run this code, I end up with these errors

Error in Error using March17Code>odefunction

Too many input arguments.

Error in March17Code>@(L,X_CH4)odefunction(L,X_CH4,totalMoleFlowRateAtBottom,k_cf,k_nf,K_c,P_initial,epsilon_CH4,fraction_H2_b,fraction_I_b,R,T,D,rho_L,sigma_L,g,nu_L) (line 78)

odefunction(L, X_CH4, totalMoleFlowRateAtBottom, k_cf, k_nf, K_c, ...

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in odearguments (line 93)

f0 = ode(t0,y0,args{:});

^^^^^^^^^^^^^^^^^^

Error in ode45 (line 104)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in March17Code (line 77)

[L, X_CH4_profile] = ode45(@(L, X_CH4) ...

^^^^^^^^^^^^^^^^^^^^^

I do not know enough MatLab to understand these errors or to troubleshoot on my own. I also have another code, Feb8Code.m in the attachments and the code is provided below. This code uses a simple Euler method that is more easy for me to understand but does not seem to be very stable. I often get nonsensical or unexpected values when I change parameters such as T, D, and P_initial under the 'Define paramters' comment.

% Define paramters

T = 1300; % Temperature [K]

D = 0.03; % Diameter [m]

P_initial = 101325 * 2; % Initial pressure (Pa)

L_points = linspace(0, 1, 1000); % Discretize L from 0 to 1

literFlowRate = 0.5; % flowrate in L/min

epsilon_CH4 = 0.8; % mole fraction of CH4 in reactor feed

epsilon_H2 = 0; % mole fraction of CH4 in reactor feed

epsilon_Ni = 0.27; % mole fraction of Ni in molten bath

V_gj = 0.1; % Drift velocity if j_g < 0.5

% Define constants

R = 8.314; % Universal gas constant (J/(mol·K))

g = 9.81; % Acceleration due to gravity (m/s^2)

A = -1.39; % Empirical constant

No = 6.02214 * 10^23; % Avogadro number

M_Ni = 58.63/1000; % Molar mass of Ni (kg/mol)

M_Bi = 208.98/1000; % Molar mass of Bi (kg/mol)

M_CH4 = 16.04/1000; % Molar mass of CH4 (kg/mol)

M_H2 = 2.02/1000; % Molar mass of H2 (kg/mol)

M_I = 39.95/1000; % Molar mass of Ar assuming Ar is I (kg/mol)

T_Melting = 1013; % Melting temperature of metal (K)

tolerance = 1e-6; % Convergence tolerance

max_iterations = 1000; % Maximum number of iterations

k_nf_0 = 14676000000; % Non-catalytic rate constant

E_nf_a = 284948; % Activation energy in J/mol

k_cf_0 = 78800; % Catalytic rate constant

E_cf_a = 208000; % Activation energy in J/mol

K = exp(13.2714 - (91204.6 / (R * T))); % Equilibrium constant

P0 = 101325; % Pressure

K_c = K * (P0 / (R * T)); % Adjust equilibrium constant

n = 1.0809;

% Initialize parameters

P = P_initial; % Initialize P

X_CH4_initial = 0; % Initialize X_CH4

X_CH4 = X_CH4_initial; % Initialize X_CH4

alphaInitialGuess = 0.5; % Initial guess for alpha_g

alpha_profile = zeros(size(L_points)); % Initialize alpha_g profile

P_profile = zeros(size(L_points)); % Initialize pressure profile

rho_g_profile = zeros(size(L_points)); % Initialize rho_g profile

Mg_profile = zeros(size(L_points)); % Initialize Mg profile for plotting

j_G_profile = zeros(size(L_points)); % Initialize j_G profile for superficial gas velocity

V_gj_plus_total_profile = zeros(size(L_points)); % Initialize array to store V_gj_plus_total values

dX_CH4_dL_profile = zeros(size(L_points));

a_g_profile = zeros(size(L_points)); % Initialize a_g profile

X_CH4_profile = zeros(size(L_points));

dvs_profile = zeros(size(L_points));

r_cat_profile = zeros(size(L_points));

r_non_profile = zeros(size(L_points));

C_CH4_profile = zeros(size(L_points));

C_H2_profile = zeros(size(L_points));

Q_G_profile = zeros(size(L_points)); % Store Q_G values for later use

% Calculated parameters

n_total = (literFlowRate / 22.4) / 60; % convert L/min to mole/s

n_CH4_b = epsilon_CH4 * n_total;

n_H2_b = epsilon_H2 * n_total;

n_I_b = (1 - epsilon_CH4 - epsilon_H2) * n_total;

totalMoleFlowRateAtBottom = n_CH4_b + n_H2_b + n_I_b;

fraction_H2_b = n_H2_b / n_CH4_b;

fraction_I_b = n_I_b / n_CH4_b;

rho_Ni = 9091 - 0.4932 * T;

rho_Bi = 10698 - 1.2064 * T;

k_cf = k_cf_0 * exp(-E_cf_a/(R * T));

k_nf = k_nf_0 * exp(-E_nf_a/(R * T));

% Density of liquid metal calculation

numerator_rho_L = epsilon_Ni * M_Ni + (1 - epsilon_Ni) * M_Bi;

denominator_rho_L = ((epsilon_Ni * M_Ni) / rho_Ni) + (((1 - epsilon_Ni) * M_Bi) / rho_Bi);

rho_L = numerator_rho_L / denominator_rho_L; % Liquid density (kg/m^3)

% Surface tension of liquid metal calculation

sigma_Bi = 0.4208 - 8.1e-5 * T;

V_Bi = M_Bi / rho_Bi;

A_Bi = 1.091 * (No^(1/3)) * (V_Bi^(2/3));

sigma_L = (sigma_Bi + ((R * T) / A_Bi) * log((1.21 * 0.57)/(1.33*0.43))); % Surface tension (N/m)

% Viscosity of liquid metal calculation

parameterB = 2.65 * (T_Melting^(1.27));

M_L = (epsilon_Ni * M_Ni) + ((1 - epsilon_Ni) * M_Bi);

numerator_parameterA = 1.7e-7 * (rho_L^(2/3)) * (T_Melting^(1/2)) * M_L^(-1/6);

denominator_parameterA = exp( parameterB / (R * T_Melting));

parameterA = numerator_parameterA / denominator_parameterA;

mew_L = parameterA * exp ( parameterB / (R * T)); % Liquid dynamic viscosity (Pa·s)

nu_L = mew_L / rho_L; % kinematic viscosity (m^2/s)

% Loop over each point in L

for idx = 1:length(L_points)

% Superficial gas velocity

j_g = ((4 * totalMoleFlowRateAtBottom * (1 + epsilon_CH4 * X_CH4)) / (pi * (D^2)))* ((R*T) / P);

% Gas density calculation

rho_g_Numerator = ((fraction_H2_b + 2 * X_CH4) * M_H2) + ((1-X_CH4) * M_CH4) + (fraction_I_b * M_I);

rho_g_Denominator = 1 + fraction_H2_b + X_CH4 + fraction_I_b;

rho_g = (rho_g_Numerator / rho_g_Denominator) * (P / (R * T));

% Gas holdup calculation

j_g_plusSubNumerator = sigma_L * g * abs(rho_L - rho_g);

j_g_plusSubDenominator = rho_L^2;

j_g_plusDenominator = j_g_plusSubNumerator / j_g_plusSubDenominator;

j_g_plus = j_g / ((j_g_plusDenominator^(0.25)));

if j_g_plus < 0.5

V_gj_plus2_SubNumerator = sigma_L * g * abs(rho_L - rho_g);

V_gj_plus2_SubDenominator = rho_L^2;

V_gj_plus2_Denominator = V_gj_plus2_SubNumerator / V_gj_plus2_SubDenominator;

V_gj_plus2 = V_gj / ((V_gj_plus2_Denominator^(0.25)));

elseif j_g_plus > 0.5

N_mewDenominator = rho_L * sigma_L * ((sigma_L/(g * abs(rho_L - rho_g)))^(0.5));

N_mew = nu_L / N_mewDenominator;

D_H_starDenominator = (sigma_L / (g * abs(rho_L - rho_g)));

D_H_star = D / D_H_starDenominator;

if N_mew <= 2.2e-3 && D_H_star <= 30

V_gj_plus2 = 0.0019 * (D_H_star*0.809) * ((rho_g/rho_L)^(-0.157)) * (N_mew^(-0.562));

elseif N_mew <= 2.2e-3 && D_H_star >= 30

V_gj_plus2 = 0.03 * ((rho_g/rho_L)^(-0.157)) * (N_mew^(-0.562));

elseif N_mew >= 2.2e-3 && D_H_star >= 30

V_gj_plus2 = 0.92 * ((rho_g/rho_L)^(-0.157));

end

end

C_0 = 1.2 - (0.2 * ((rho_g / rho_L)^(0.5)));

alpha = alphaInitialGuess; % Initial guess for alpha_g

for iteration = 1:max_iterations

V_gj_plus1 = (2^(0.5)) * ((1- alpha)^1.75);

V_gj_plus = (V_gj_plus1 * exp(A * j_g_plus)) + (V_gj_plus2 * (1 - exp(A * j_g_plus)));

alpha_new = j_g_plus / ((C_0 * j_g_plus) + V_gj_plus);

% Check for convergence

if abs(alpha_new - alpha) < tolerance

alpha_profile(idx) = alpha_new;

break;

end

% Update alpha_g for next iteration

alpha = alpha_new;

% Stop if maximum iterations are reached

if iteration == max_iterations

warning('Maximum iterations reached at L = . Alpha_g may not have converged.');

end

end

% Specific interfacial area calculation

N_Bo = (g * (D^2) * rho_L) / sigma_L;

N_Ga = (g * (D^3)) / (nu_L^2);

N_Fr = j_g / ((g * D)^(0.5));

d_vs = D * 26 * (N_Bo^(-0.5)) * (N_Ga^(-0.12)) * (N_Fr^(-0.12));

a_g = 6 / d_vs;

a_g_profile(idx) = a_g;

% Update pressure at current L

dP_dL = -((rho_L * (1 - alpha)) + (rho_g * alpha)) * g;

P = P + dP_dL * (L_points(2) - L_points(1)); % Incremental pressure update

P_profile(idx) = P; % Store X_CH4 value for the current L

% X_CH4 as a function of L

term0 = (alpha * pi * D^2) / 4;

term1_Numerator = a_g * k_cf * (1 - X_CH4);

term1_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4));

term1 = (term1_Numerator / term1_Denominator) * (P / (R * T));

subTerm2_Numerator = 1 - X_CH4;

subTerm2_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4));

subTerm2 = (subTerm2_Numerator / subTerm2_Denominator) * (P / (R * T));

term2 = k_nf * (subTerm2^n);

subTerm3_Numerator = n_CH4_b * ((fraction_H2_b + 2 * X_CH4)^2);

subTerm3_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4)) * (1 - X_CH4) * K_c;

subTerm3 = (subTerm3_Numerator / subTerm3_Denominator) * (P / (R * T));

term3 = 1 - subTerm3;

dX_CH4_dL = term0*(term1 + term2)*(term3);

X_CH4 = X_CH4 + dX_CH4_dL * (L_points(2) - L_points(1));

X_CH4_profile(idx) = X_CH4; % Store X_CH4 value for the current L

end

figure;

plot(L_points, X_CH4_profile, 'b-', 'LineWidth', 2);

xlabel('Reactor Length (L) [m]');

ylabel('CH_4 Conversion (X_{CH4})');

title('CH_4 Conversion Profile along Reactor Length');

grid on;

legend('X_{CH4} vs L');

figure;

plot(L_points, P_profile, 'b-', 'LineWidth', 2);

xlabel('Reactor Length (L) [m]');

ylabel('Pressure (Pa)');

title('P Profile along Reactor Length');

grid on;

legend('P vs L');

Any help or insight is much appreciated.

Answer 1

[Sam Chak]

Hi @Tiek Aun

I am unfamiliar with the dynamics of the pyrolysis reactor; however, I can point out two mistakes in the code.

First, if the reactor has many parameters that are treated as constants, it may be advisable to declare them in the odefunction(). Otherwise, the input argument will consist of a long chain of parameters. The error message indicates that the input argument of odefunction() is missing the specified variables. In fact, ten parameters are missing:

M_H2, M_CH4, M_I, V_gj, alphaInitialGuess, max_iterations, A, tolerance, n, n_CH4_b.

The second mistake is that there is another ODE, , embedded within the ODE for . You attempt to solve using the Euler method, which is prone to the accumulation of errors at each L step as the method progresses. It is possible to define as the second state equation for ode45() to solve.

The following code now runs, but I am unsure if everything else is correct.

% Define paramters
T = 1150; % Temperature [K]
D = 0.03; % Diameter [m]
P_initial = 101325 * 3; % Initial pressure (Pa)
L_span = [0 2]; % Discretize L from 0 to 1
literFlowRate = 0.5; % flowrate in L/min
epsilon_CH4 = 0.8; % mole fraction of CH4 in reactor feed
epsilon_H2 = 0; % mole fraction of CH4 in reactor feed
epsilon_Ni = 0.27; % mole fraction of Ni in molten bath
V_gj = 0.1; % Drift velocity if j_g < 0.5

% Define constants
R = 8.314; % Universal gas constant (J/(mol·K))
g = 9.81; % Acceleration due to gravity (m/s^2)
A = -1.39; % Empirical constant
No = 6.02214 * 10^23; % Avogadro number
M_Ni = 58.63/1000; % Molar mass of Ni (kg/mol)
M_Bi = 208.98/1000; % Molar mass of Bi (kg/mol)
M_CH4 = 16.04/1000; % Molar mass of CH4 (kg/mol)
M_H2 = 2.02/1000; % Molar mass of H2 (kg/mol)
M_I = 39.95/1000; % Molar mass of Ar assuming Ar is I (kg/mol)
T_Melting = 1013; % Melting temperature of metal (K)
tolerance = 1e-6; % Convergence tolerance
max_iterations = 1000; % Maximum number of iterations
k_nf_0 = 14676000000; % Non-catalytic rate constant
E_nf_a = 284948; % Activation energy in J/mol
k_cf_0 = 78800; % Catalytic rate constant
E_cf_a = 208000; % Activation energy in J/mol
K = exp(13.2714 - (91204.6 / (R * T))); % Equilibrium constant
P0 = 101325; % Pressure
K_c = K * (P0 / (R * T)); % Adjust equilibrium constant
n = 1.0809;

% Initialize parameters
P = P_initial; % Initialize P
X_CH4_initial = 0; % Initialize X_CH4
alphaInitialGuess = 0.5; % Initial guess for alpha_g

% Calculated parameters
n_total = (literFlowRate / 22.4) / 60; % convert L/min to mole/s
n_CH4_b = epsilon_CH4 * n_total;
n_H2_b = epsilon_H2 * n_total;
n_I_b = (1 - epsilon_CH4 - epsilon_H2) * n_total;
totalMoleFlowRateAtBottom = n_CH4_b + n_H2_b + n_I_b;
fraction_H2_b = n_H2_b / n_CH4_b;
fraction_I_b = n_I_b / n_CH4_b;
rho_Ni = 9091 - 0.4932 * T;
rho_Bi = 10698 - 1.2064 * T;
k_cf = k_cf_0 * exp(-E_cf_a/(R * T));
k_nf = k_nf_0 * exp(-E_nf_a/(R * T));

% Density of liquid metal calculation
numerator_rho_L = epsilon_Ni * M_Ni + (1 - epsilon_Ni) * M_Bi;
denominator_rho_L = ((epsilon_Ni * M_Ni) / rho_Ni) + (((1 - epsilon_Ni) * M_Bi) / rho_Bi);
rho_L = numerator_rho_L / denominator_rho_L; % Liquid density (kg/m^3)
% Surface tension of liquid metal calculation
sigma_Bi = 0.4208 - 8.1e-5 * T;
V_Bi = M_Bi / rho_Bi;
A_Bi = 1.091 * (No^(1/3)) * (V_Bi^(2/3));
sigma_L = (sigma_Bi + ((R * T) / A_Bi) * log((1.21 * 0.57)/(1.33*0.43))); % Surface tension (N/m)

% Viscosity of liquid metal calculation
parameterB = 2.65 * (T_Melting^(1.27));
M_L = (epsilon_Ni * M_Ni) + ((1 - epsilon_Ni) * M_Bi);
numerator_parameterA = 1.7e-7 * (rho_L^(2/3)) * (T_Melting^(1/2)) * M_L^(-1/6);
denominator_parameterA = exp( parameterB / (R * T_Melting));
parameterA = numerator_parameterA / denominator_parameterA;
mew_L = parameterA * exp ( parameterB / (R * T)); % Liquid dynamic viscosity (Pa·s)
nu_L = mew_L / rho_L; % kinematic viscosity (m^2/s)

[L, X_CH4_profile] = ode45(@(L, state) odefunction(L, state, totalMoleFlowRateAtBottom, k_cf, k_nf, K_c, epsilon_CH4, fraction_H2_b, fraction_I_b, R, T, D, rho_L, sigma_L, g, nu_L, M_H2, M_CH4, M_I, V_gj, alphaInitialGuess, max_iterations, A, tolerance, n, n_CH4_b), L_span, [X_CH4_initial; P_initial]);

% Loop over each point in L
function dX_CH4_dL = odefunction(L, state, totalMoleFlowRateAtBottom, k_cf, k_nf, K_c, epsilon_CH4, fraction_H2_b, fraction_I_b, R, T, D, rho_L, sigma_L, g, nu_L, M_H2, M_CH4, M_I, V_gj, alphaInitialGuess, max_iterations, A, tolerance, n, n_CH4_b)

    X_CH4   = state(1);
    P       = state(2);
   
    % Superficial gas velocity
    j_g = ((4 * totalMoleFlowRateAtBottom * (1 + epsilon_CH4 * X_CH4)) / (pi * (D^2)))* ((R*T) / P);
    
    % Gas density calculation
    rho_g_Numerator = ((fraction_H2_b + 2 * X_CH4) * M_H2) + ((1-X_CH4) * M_CH4) + (fraction_I_b * M_I);
    rho_g_Denominator = 1 + fraction_H2_b + X_CH4 + fraction_I_b;
    rho_g = (rho_g_Numerator / rho_g_Denominator) * (P / (R * T));
    
    % Gas holdup calculation
    j_g_plusSubNumerator = sigma_L * g * abs(rho_L - rho_g);
    j_g_plusSubDenominator = rho_L^2;
    j_g_plusDenominator = j_g_plusSubNumerator /  j_g_plusSubDenominator;
    j_g_plus = j_g / ((j_g_plusDenominator^(0.25)));
    
    if j_g_plus < 0.5
        V_gj_plus2_SubNumerator = sigma_L * g * abs(rho_L - rho_g);
        V_gj_plus2_SubDenominator = rho_L^2;
        V_gj_plus2_Denominator = V_gj_plus2_SubNumerator /  V_gj_plus2_SubDenominator;
        V_gj_plus2 = V_gj / ((V_gj_plus2_Denominator^(0.25)));
    elseif j_g_plus > 0.5
        N_mewDenominator = rho_L * sigma_L * ((sigma_L/(g * abs(rho_L - rho_g)))^(0.5));
        N_mew = nu_L / N_mewDenominator;
        D_H_starDenominator = (sigma_L / (g * abs(rho_L - rho_g)));
        D_H_star = D / D_H_starDenominator;
        if N_mew <= 2.2e-3 && D_H_star <= 30
            V_gj_plus2 = 0.0019 * (D_H_star*0.809) * ((rho_g/rho_L)^(-0.157)) * (N_mew^(-0.562));
        elseif N_mew <= 2.2e-3 && D_H_star >= 30
            V_gj_plus2 = 0.03 * ((rho_g/rho_L)^(-0.157)) * (N_mew^(-0.562));
        elseif N_mew >= 2.2e-3 && D_H_star >= 30
            V_gj_plus2 = 0.92 * ((rho_g/rho_L)^(-0.157));
        end
    end
    
    C_0 = 1.2 - (0.2 * ((rho_g / rho_L)^(0.5)));
    
    alpha = alphaInitialGuess; % Initial guess for alpha_g
    
    for iteration = 1:max_iterations
        V_gj_plus1 = (2^(0.5)) * ((1- alpha)^1.75);
        V_gj_plus = (V_gj_plus1 * exp(A * j_g_plus)) + (V_gj_plus2 * (1 - exp(A * j_g_plus)));
        alpha_new = j_g_plus / ((C_0 * j_g_plus) + V_gj_plus);
         % Check for convergence
            if abs(alpha_new - alpha) < tolerance
              
                break;
            end
            % Update alpha_g for next iteration
            alpha = alpha_new;
            % Stop if maximum iterations are reached
            if iteration == max_iterations
                warning('Maximum iterations reached at L = . Alpha_g may not have converged.');
            end
    end
    
    % Specific interfacial area calculation
    N_Bo = (g * (D^2) * rho_L) / sigma_L;
    N_Ga = (g * (D^3)) / (nu_L^2);
    N_Fr = j_g / ((g * D)^(0.5));
    
    d_vs = D * 26 * (N_Bo^(-0.5)) * (N_Ga^(-0.12)) * (N_Fr^(-0.12));
    a_g = 6 / d_vs;
    
    % Update pressure at current L
    dP_dL = -((rho_L * (1 - alpha)) + (rho_g * alpha)) * g;
    % P     = P + dP_dL * (L(2) - L(1));  % Incremental pressure update
 
    % X_CH4 as a function of L 
    term0 = (alpha * pi * D^2) / 4;
    
    term1_Numerator = a_g * k_cf * (1 - X_CH4);
    term1_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4));
    term1 = (term1_Numerator / term1_Denominator) * (P / (R * T));
    
    subTerm2_Numerator = 1 - X_CH4;
    subTerm2_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4));
    subTerm2 = (subTerm2_Numerator / subTerm2_Denominator) * (P / (R * T));
    term2 = k_nf * (subTerm2^n);
    
    subTerm3_Numerator = n_CH4_b * ((fraction_H2_b + 2 * X_CH4)^2);
    subTerm3_Denominator = totalMoleFlowRateAtBottom * (1 + (epsilon_CH4 * X_CH4)) * (1 - X_CH4) * K_c;
    subTerm3 = (subTerm3_Numerator / subTerm3_Denominator) * (P / (R * T));
    term3 = 1 - subTerm3;
    
    dX_CH4_dL(1) = term0*(term1 + term2)*(term3);
    dX_CH4_dL(2) = dP_dL;
    dX_CH4_dL    = [dX_CH4_dL(1)
                    dX_CH4_dL(2)];
    
end

% Plot CH4 Conversion along Reactor Length
figure;
plot(L, X_CH4_profile(:,1), 'b-', 'LineWidth', 2);
xlabel('Reactor Length (L) [m]');
ylabel('CH_4 Conversion (X_{CH4})');
title('CH_4 Conversion Profile along Reactor Length');
grid on;
legend('X_{CH4} vs L');

[Sam Chak]

The first error message, 'Too many input arguments,' refers to the error in odefunction(). Upon inspection, you will find that P_initial is missing. Since it serves as the initial value for ​ inside odefunction(), every time the ode45 solver runs, it will repeatedly overwrite the evolved P state in with P_initial.

[Tiek Aun]

Thank you for the answer and comment. I think I understand better now. Much appreciated.