Not enough input arguments.

Dear everyone,

I hope you are having a great day.

I found a code attached to a research whom his developer advices that is working regarding modelling biochemical process. However. when trying running it, I found some problems.

One is when I run the following function

function y1 = ADDE(t,y)
y1 = zeros(size(y));

and so on giving equations for y1(elements)

comes an error of "Not enough input arguments".

So should I define y first? because it is not performed in the orginial code.

Thanks

[Star Strider]

That function appears to be designed to use the MATLAB numeric differential equation integrators (such as ode45, or with a biochemical model, more likely ode15s).

If you use it in that context, it should work correctly, without any changes being made to it.

Check the paper to be certain that is how it is intended to be used.

.

[Mohamed Eisa]

Yes sir, it's exactly as you are saying (ode15s). Without changing anything and by sticking completely to all the guidelines of the paper, it shows the error of not enough input arguments. It's the first time I am being objected to such a code, therefore I am a little bit struggling.

[Mohamed Eisa]

The model consists of 3 parts, first part is code 1(first m file) to take the general input variables values from excel, second part (second m file) which contains all the differential equations and calculations and which I get this error from it while running it, third part (third m file) contains the ode15s function and plotting and it gives also the same problem.

I couldn't find any help so far regarding it and I am trying so hard to configure what the reason of the error might be so any help of you is totally appreciated! Many thanks for your time.

[Star Strider]

My pleasure.

It would be helpful if you could post / attach a .pdf version of the complete paper (including any necessary supplements, since they usually contain the appendices with data and code).

[Mohamed Eisa]

Dear Gentleman,

Firstly, I can't thank you enough for your consideration as I am stuck on that with no one around experienced enough to ask. I am attaching the codes files that I copied and the excel sheet of the inputs which took me around 2 hrs to copy from the pdf.

Also, here is a driect link of the Thesis pdf ( https://eprints.soton.ac.uk/375082/2/HHN_Thesis_FINAL_Feb_2017_rechecked.pdf ) as I could not upload it for exceeding the limit of 5MB.

Thanks a lot, I really realy appreciate it.

Yours,

Mohamed

[Walter Roberson]

To use that code you first have to run globalvariables in order to read data files and set up the global variables.

But after that, you need code that runs the ode15s function, and you have not posted that code.

Do not run ADDE by itself. You would need a call that looked more like

ode15s(@ADDE, tspan,initial_conditions)

and notice the @ there -- it is important! If you were to run

ode15s(ADDE, tspan,initial_conditions)

without the @ then you would get the error message you show.

[Mohamed Eisa]

Dear Mr. Roberson,

Many thanks for your time and advice.

I usually run first the global variables as suggested and then there are two codes (ADDE with the calculations itself and AD with ode15s equation); anytime I try to run each one of those; I got the same error (not enough inputs).

Now, I tried @ before ADDE in the ode15s function and it is the same error.

Thanks for your time!

Yours,

Mohamed

[Walter Roberson]

Please post the code you use to make the ode15s call.

[Mohamed Eisa]

Thanks!

That is for global variables

format long
clear all

global q
global V_dig
global V_liq
global V_gas
global tspan
global u
global maxx

global S_suf
global S_aaf
global S_faf
global S_vaf
global S_buf
global S_prof
global S_acf
global S_h2f
global S_ch4f
global S_ICf
global S_INf
global S_If
global X_cf
global X_chf
global X_prf
global X_lif
global X_suf
global X_aaf
global X_faf
global X_c4f 
global X_prof
global X_acf
global X_ac2f
global X_h2f
global X_If
global S_cat_ionf
global S_an_ionf
global S_va_ionf
global S_bu_ionf
global S_pro_ionf
global S_ac_ionf
global S_hco3_ionf
global S_nh3f
global S_gas_h2f
global S_gas_ch4f
global S_gas_co2f
global S_h_ionf
global f_sI_xc
global f_xI_xc
global f_ch_xc
global f_pr_xc
global f_li_xc
global f_fa_li
global f_h2_su
global f_bu_su
global f_pro_su
global f_ac_su
global f_h2_aa
global f_va_aa
global f_bu_aa
global f_pro_aa
global f_ac_aa
global C_xc
global C_sI

global Y_su
global Y_aa
global Y_fa
global Y_c4
global Y_pro
global Y_ac
global Y_ac2
global Y_h2 
global k_dis
global k_hyd_ch
global k_hyd_pr
global k_hyd_li
global k_m_su
global k_m_aa
global k_m_fa
global k_m_c4
global k_m_pro
global k_m_ac
global k_m_ac2
global k_m_h2
global k_dec_Xsu
global k_dec_Xaa
global k_dec_Xfa
global k_dec_Xc4
global k_dec_Xpro
global k_dec_Xac
global k_dec_Xac2
global k_dec_Xh2
global K_S_IN
global K_S_su
global K_S_aa
global K_S_fa
global K_Ih2_fa
global K_S_pro
global K_Ih2_pro
global K_S_ac
global K_S_ac2
global K_I_nh3
global K_S_c4
global K_S_h2
global K_Ih2_c4
global kLa
global K_H_h2o_base
global K_H_co2_base
global K_H_ch4_base
global K_H_h2_base
global k_P
global P_atm
global T_base
global T_op
global R
global pK_w_base
global pK_a_va_base
global pK_a_bu_base
global pK_a_pro_base
global pK_a_ac_base
global pK_a_co2_base
global pK_a_IN_base
global k_A_Bva
global k_A_Bbu
global k_A_Bpro
global k_A_Bac
global k_A_Bco2
global k_A_BIN
global pH_UL_h2
global pH_LL_h2
global pH_UL_aa
global pH_LL_aa
global pH_UL_ac
global pH_LL_ac
global pH_UL_ac2
global pH_LL_ac2
global K_Ih2_ac
global pHLim_aa
global pHLim_ac
global pHLim_ac2
global pHLim_h2
global k_aa
global k_ac
global k_ac2
global k_h2
global I11a;
global I11b;
global I18a;
global I18b;
u = xlsread('ADM1data','Digesterconfig','K3:K39'); % Initial conditions
Inputs = xlsread('ADM1data','Digesterconfig','F3:F39');
Digesterconfig = xlsread('ADM1data','Digesterconfig','B3:B9');
Fraction = xlsread('ADM1data','Biostoic','C3 : C17');
Carbonstoichiometries = xlsread('ADM1data','Biostoic','C21 : C35');
Nitrogenstoichiometries = xlsread('ADM1data','Biostoic','C39 : C42');
Yielduptakecomponents = xlsread('ADM1data','Biochemrate','C3:C10');
Dishydcoefficients = xlsread('ADM1data','Biochemrate','C13:C32');
Halfsaturatecoefficients = xlsread('ADM1data','Biochemrate','C34:C47');
Acidgasparameters = xlsread('ADM1data','Biochemrate','C50:C80');
q = Digesterconfig(1);
V_dig = Digesterconfig(2);
V_liq = Digesterconfig(3);
V_gas = Digesterconfig(4);
tspan = [Digesterconfig(6) Digesterconfig(7)];
maxx = Digesterconfig(7);
S_suf = Inputs(1);
S_aaf = Inputs(2);
S_faf = Inputs(3);
S_vaf = Inputs(4);
S_buf = Inputs(5);
S_prof = Inputs(6);
S_acf = Inputs(7);
S_h2f = Inputs(8);
S_ch4f = Inputs(9);
S_ICf = Inputs(10);
S_INf = Inputs(11);
S_If = Inputs(12);
X_cf = Inputs(13);
X_chf = Inputs(14);
X_prf = Inputs(15);
X_lif = Inputs(16);
X_suf = Inputs(17);
X_aaf = Inputs(18);
X_faf = Inputs(19);
X_c4f = Inputs(20);
X_prof = Inputs(21);
X_acf = Inputs(22);
X_h2f = Inputs(23);
X_If = Inputs(24);
S_cat_ionf = Inputs(25);
S_an_ionf = Inputs(26);
S_va_ionf = Inputs(27);
S_bu_ionf = Inputs(28);
S_pro_ionf = Inputs(29);
S_ac_ionf = Inputs(30);
S_hco3_ionf = Inputs(31);
S_nh3f = Inputs(32);
S_gas_h2f = Inputs(33);
S_gas_ch4f = Inputs(34);
S_gas_co2f = Inputs(35);
S_h_ionf = Inputs(36);
X_ac2f = Inputs(37);
f_sI_xc = Fraction (1,1);
f_xI_xc = Fraction (2,1);
f_ch_xc = Fraction (3,1);
f_pr_xc = Fraction (4,1);
f_li_xc = Fraction (5,1);
f_fa_li = Fraction (6,1);
f_h2_su = Fraction (7,1);
f_bu_su = Fraction (8,1);
f_pro_su = Fraction (9,1);
f_ac_su = Fraction (10,1);
f_h2_aa = Fraction (11,1);
f_va_aa = Fraction (12,1);
f_bu_aa = Fraction (13,1);
f_pro_aa = Fraction (14,1);
f_ac_aa = Fraction (15,1);
C_xc = Carbonstoichiometries (1);
C_sI = Carbonstoichiometries (2);
C_ch = Carbonstoichiometries (3);
C_pr = Carbonstoichiometries (4);
C_li = Carbonstoichiometries (5);
C_xI = Carbonstoichiometries (6);
C_su = Carbonstoichiometries (7);
C_aa = Carbonstoichiometries (8);
C_fa = Carbonstoichiometries (9);
C_bu = Carbonstoichiometries (10);
C_pro = Carbonstoichiometries (11);
C_ac = Carbonstoichiometries (12);
C_bac = Carbonstoichiometries (13);
C_va = Carbonstoichiometries (14);
C_ch4 = Carbonstoichiometries (15);
N_xc = Nitrogenstoichiometries (1);
N_I = Nitrogenstoichiometries (2);
N_aa = Nitrogenstoichiometries (3);
N_bac = Nitrogenstoichiometries (4);
Y_su = Yielduptakecomponents (1);
Y_aa = Yielduptakecomponents (2);
Y_fa = Yielduptakecomponents (3);
Y_c4 = Yielduptakecomponents (4);
Y_pro = Yielduptakecomponents (5);
Y_ac = Yielduptakecomponents (6);
Y_h2 = Yielduptakecomponents (7);
Y_ac2 = Yielduptakecomponents (8);
k_dis = Dishydcoefficients (1);
k_hyd_ch = Dishydcoefficients (2);
k_hyd_pr = Dishydcoefficients (3);
k_hyd_li = Dishydcoefficients (4);
k_m_su = Dishydcoefficients (5);
k_m_aa = Dishydcoefficients (6);
k_m_fa = Dishydcoefficients (7);
k_m_c4 = Dishydcoefficients (8);
k_m_pro = Dishydcoefficients (9);
k_m_ac = Dishydcoefficients (10);
k_m_h2 = Dishydcoefficients (11);
k_dec_Xsu = Dishydcoefficients (12);
k_dec_Xaa = Dishydcoefficients (13);
k_dec_Xfa = Dishydcoefficients (14);
k_dec_Xc4 = Dishydcoefficients (15);
k_dec_Xpro = Dishydcoefficients (16);
k_dec_Xac = Dishydcoefficients (17);
k_dec_Xh2 = Dishydcoefficients (18);
k_m_ac2 = Dishydcoefficients (19);
k_dec_Xac2 = Dishydcoefficients (20);
K_S_IN = Halfsaturatecoefficients (1);
K_S_su = Halfsaturatecoefficients (2);
K_S_aa = Halfsaturatecoefficients (3);
K_S_fa = Halfsaturatecoefficients (4);
K_Ih2_fa = Halfsaturatecoefficients (5);
K_S_pro = Halfsaturatecoefficients (6);
K_Ih2_pro = Halfsaturatecoefficients (7);
K_S_ac = Halfsaturatecoefficients (8);
K_I_nh3 = Halfsaturatecoefficients (9);
K_S_c4 = Halfsaturatecoefficients (10);
K_S_h2 = Halfsaturatecoefficients (11);
K_Ih2_c4 = Halfsaturatecoefficients (12);
K_S_ac2 = Halfsaturatecoefficients (13);
K_Ih2_ac = Halfsaturatecoefficients (14);
kLa = Acidgasparameters (1);
K_H_h2o_base = Acidgasparameters (2);
K_H_co2_base = Acidgasparameters (3);
K_H_ch4_base = Acidgasparameters (4);
K_H_h2_base = Acidgasparameters (5);
k_P = Acidgasparameters (6);
P_atm = Acidgasparameters (7);
T_base = Acidgasparameters (8);
T_op = Acidgasparameters (9);
R = Acidgasparameters (10);
pK_w_base = Acidgasparameters (11);
pK_a_va_base = Acidgasparameters (12);
pK_a_bu_base = Acidgasparameters (13);
pK_a_pro_base = Acidgasparameters (14);
pK_a_ac_base = Acidgasparameters (15);
pK_a_co2_base = Acidgasparameters (16);
pK_a_IN_base = Acidgasparameters (17);
k_A_Bva = Acidgasparameters (18);
k_A_Bbu = Acidgasparameters (19);
k_A_Bpro = Acidgasparameters (20);
k_A_Bac = Acidgasparameters (21);
k_A_Bco2 = Acidgasparameters (22);
k_A_BIN = Acidgasparameters (23);
pH_UL_h2 = Acidgasparameters (24);
pH_LL_h2 = Acidgasparameters (25);
pH_UL_aa = Acidgasparameters (26);
pH_LL_aa = Acidgasparameters (27);
pH_UL_ac = Acidgasparameters (28);
pH_LL_ac = Acidgasparameters (29);
pH_UL_ac2 = Acidgasparameters (30);
pH_LL_ac2 = Acidgasparameters (31);
pHLim_aa = 10^(-(pH_UL_aa + pH_LL_aa)/2.0);
pHLim_ac = 10^(-(pH_UL_ac + pH_LL_ac)/2.0);
pHLim_ac2 = 10^(-(pH_UL_ac2 + pH_LL_ac2)/2.0);
pHLim_h2 = 10^(-(pH_UL_h2 + pH_LL_h2)/2.0);
k_aa = 3.0/(pH_UL_aa-pH_LL_aa);
k_ac = 3.0/(pH_UL_ac-pH_LL_ac);
k_ac2 = 3.0/(pH_UL_ac2-pH_LL_ac2);
k_h2 = 3.0/(pH_UL_h2-pH_LL_h2);
I11a = 1;
I11b = 1;
I18a = 1;
I18b = 1;
clear Inputs;
clear Digesterconfig;
clear Fraction;
clear Carbonstoichiometries;
clear Nitrogenstoichiometries;
clear Yielduptakecomponents;
clear Dishydcoefficients;
clear Halfsaturatecoefficients;
clear Acidgasparameters;
save saveddata

That is for all the calculations and ODE equations (ADDE.m)

function y1 = ADDE(t,y)
y1 = zeros(size(y));

format long
global q
global V_liq % Volume of liquid part
global V_gas % Volume of gas space
global S_suf
global S_aaf
global S_faf
global S_vaf
global S_buf
global S_prof
global S_acf
global S_h2f
global S_ch4f
global S_ICf
global S_INf
global S_If
global X_cf
global X_chf
global X_prf
global X_lif
global X_suf
global X_aaf
global X_faf
global X_c4f
global X_prof
global X_acf
global X_ac2f
global X_h2f
global X_If
global S_cat_ionf
global S_an_ionf
global f_sI_xc
global f_xI_xc
global f_ch_xc
global f_pr_xc
global f_li_xc
global f_fa_li
global f_h2_su
global f_bu_su
global f_pro_su
global f_ac_su
global f_h2_aa
global f_va_aa
global f_bu_aa
global f_pro_aa
global f_ac_aa
global C_xc
global C_sI
global C_ch
global C_pr
global C_li
global C_xI
global C_su
global C_aa
global C_fa
global C_bu
global C_pro
global C_ac
global C_bac
global C_va
global C_ch4
global N_xc
global N_I
global N_aa
global N_bac
global Y_su
global Y_aa
global Y_fa
global Y_c4
global Y_pro
global Y_ac
global Y_ac2
global Y_h2
global k_dis
global k_hyd_ch
global k_hyd_pr
global k_hyd_li
global k_m_su
global k_m_aa
global k_m_fa
global k_m_c4
global k_m_pro
global k_m_ac
global k_m_ac2
global k_m_h2
global k_dec_Xsu
global k_dec_Xaa
global k_dec_Xfa
global k_dec_Xc4
global k_dec_Xpro
global k_dec_Xac
global k_dec_Xac2
global k_dec_Xh2
global K_S_IN
global K_S_su
global K_S_aa
global K_S_fa
global K_Ih2_fa
global K_S_pro
global K_Ih2_pro
global K_S_ac
global K_S_ac2
global K_I_nh3
global K_S_c4
global K_S_h2
global K_Ih2_c4
global K_Ih2_ac
global kLa
global K_H_h2o_base
global K_H_co2_base
global K_H_ch4_base
global K_H_h2_base
global k_P
global P_atm
global T_base
global T_op
global R
global pK_w_base
global pK_a_va_base
global pK_a_bu_base
global pK_a_pro_base
global pK_a_ac_base
global pK_a_co2_base
global pK_a_IN_base
global k_A_Bva
global k_A_Bbu
global k_A_Bpro
global k_A_Bac
global k_A_Bco2
global k_A_BIN
global pHLim_aa;
global pHLim_ac;
global pHLim_ac2;
global pHLim_h2;
global k_aa;
global k_ac;
global k_ac2;
global k_h2;
global I11a;
global I11b;
global I18a;
global I18b;
global inhib11b;
global inhib56;
global inhib7;
global inhib89;
global inhib10;
global inhib11;
global inhib12;
global Itec4;
global Itepro;
global K_a_va
global K_a_bu
global K_a_pro
global K_a_ac
global K_a_co2
global K_a_IN
global K_w
global K_H_h2
global K_H_ch4
global K_H_co2
global p_gas_h2o
global factor
global P_gas
global p_gas_h2
global p_gas_ch4
global p_gas_co2
global q_gas
factor = (1.0/T_base - 1.0/T_op)/(100*R);
K_a_va =10^(-pK_a_va_base);
K_a_bu = 10^(-pK_a_bu_base);
K_a_pro = 10^(-pK_a_pro_base);
K_a_ac = 10^(-pK_a_ac_base);
K_a_co2 = (10^(-pK_a_co2_base))*exp(7646.0*factor);
K_a_IN = 10^(-pK_a_IN_base)*exp(51965.0*factor); 
K_w = 10^(-pK_w_base)*exp(55900.0*factor); % T adjustment for K_w
K_H_h2 = K_H_h2_base*exp(-4180.0*factor); %/* T adjustment for K_H_h2
K_H_ch4 = K_H_ch4_base*exp(-14240.0*factor);
K_H_co2 = K_H_co2_base*exp(-19410.0*factor);
r_A_4 = ( k_A_Bva*(y(27)*(K_a_va + y(36)) - K_a_va*y(4)) );
r_A_5 = ( k_A_Bbu*(y(28)*(K_a_bu + y(36)) - K_a_bu*y(5)) );
r_A_6 = ( k_A_Bpro*(y(29)*(K_a_pro + y(36))- K_a_pro*y(6)) );
r_A_7 = ( k_A_Bac*(y(30)*(K_a_ac + y(36)) - K_a_ac*y(7)) );
r_A_10 =( k_A_Bco2*(y(31)*(K_a_co2 + y(36)) - K_a_co2*y(10)) ); %
r_A_11 =( k_A_BIN*(y(32)*(K_a_IN + y(36)) - K_a_IN*y(11)) ); %
r_T_8 = ( kLa*(y(8)-16*K_H_h2*( y(33)*R*T_op/16.0 )) );
r_T_9 = ( kLa*(y(9)-64*K_H_ch4*( y(34)*R*T_op/64.0 )) );
r_T_10 =( kLa*(y(10)-y(31)-K_H_co2*( y(35)*R*T_op )) );
stoich1 = -C_xc+f_sI_xc*C_sI+f_ch_xc*C_ch+f_pr_xc*C_pr+f_li_xc*C_li+f_xI_xc*C_xI;
stoich2 = -C_ch+C_su;
stoich3 = -C_pr+C_aa;
stoich4 = -C_li+(1.0-f_fa_li)*C_su+f_fa_li*C_fa;
stoich5 = -C_su+(1.0-Y_su)*(f_bu_su*C_bu+f_pro_su*C_pro+f_ac_su*C_ac)+Y_su*C_bac;
stoich6 = -C_aa+(1.0-Y_aa)*(f_va_aa*C_va+f_bu_aa*C_bu+f_pro_aa*C_pro+f_ac_aa*C_ac)+Y_aa*C_bac;
stoich7 = -C_fa+(1.0-Y_fa)*0.7*C_ac+Y_fa*C_bac;
stoich8 = -C_va+(1.0-Y_c4)*0.54*C_pro+(1.0-Y_c4)*0.31*C_ac+Y_c4*C_bac;
stoich9 = -C_bu+(1.0-Y_c4)*0.8*C_ac+Y_c4*C_bac;
stoich10 = -C_pro+(1.0-Y_pro)*0.57*C_ac+Y_pro*C_bac;
stoich11 = -C_ac+(1.0-Y_ac)*C_ch4+Y_ac*C_bac;
stoich11b = -C_ac+Y_ac2*C_bac;
stoich12 = (1.0-Y_h2)*C_ch4+Y_h2*C_bac;
stoich13 = -C_bac+C_xc;
I_pH_aa = pHLim_aa^k_aa /(y(36)^k_aa + pHLim_aa^k_aa);
I_pH_ac = pHLim_ac^k_ac /(y(36)^k_ac + pHLim_ac^k_ac);
I_pH_ac2 = pHLim_ac2^k_ac2 /(y(36)^k_ac2 + pHLim_ac2^k_ac2);
I_pH_h2 = pHLim_h2^k_h2 /(y(36)^k_h2 + pHLim_h2^k_h2);
I_IN_lim = 1.0/(1.0+K_S_IN/y(11));
I_h2_fa = 1.0/(1.0+y(8)/K_Ih2_fa);
I_h2_c4 = 1.0/(1.0+y(8)/K_Ih2_c4);
I_h2_pro = 1.0/(1.0+y(8)/K_Ih2_pro);
I_h2_ac = 1.0/(1.0+y(8)/K_Ih2_ac);
I_nh3 = 1.0/(1.0+y(32)/K_I_nh3);
Itec4 = 1;
Itepro = 1;
inhib56 = I_pH_aa*I_IN_lim;
inhib7 = inhib56*I_h2_fa;
inhib89 = inhib56*I_h2_c4*Itec4;
inhib10 = inhib56*I_h2_pro*Itepro;
inhib11 = I_pH_ac*I_IN_lim*I_nh3*I11a;
inhib11b = I_pH_ac2*I_IN_lim*I_h2_ac*I11b;
inhib12 = I_pH_h2*I_IN_lim;
ro1 = k_dis*y(13);
ro2 = k_hyd_ch*y(14);
ro3 = k_hyd_pr*y(15);
ro4 = k_hyd_li*y(16);
ro5 = k_m_su*(y(1)/(y(1)+K_S_su))*y(17)*inhib56;
ro6 = k_m_aa*(y(2)/(K_S_aa+y(2)))*y(18)*inhib56;
ro7 = k_m_fa*(y(3)/(K_S_fa+y(3)))*y(19)*inhib7;
ro8 = k_m_c4*(y(4)/(K_S_c4+y(4)))*y(20)*(y(4)/(y(5)+y(4)+1e-6))*inhib89;
ro9 = k_m_c4*(y(5)/(K_S_c4+y(5)))*y(20)*(y(5)/(y(4)+y(5)+1e-6))*inhib89;
ro10 = k_m_pro*(y(6)/(K_S_pro+y(6)))*y(21)*inhib10;
ro11 = k_m_ac*(y(7)/(K_S_ac+y(7)))*y(22)*inhib11;
ro11b = k_m_ac2*(y(7)/(K_S_ac2+y(7)))*y(37)*inhib11b; % Acetate Oxidation
ro12 = k_m_h2*(y(8)/(K_S_h2+y(8)))*y(23)*inhib12;
ro13 = k_dec_Xsu*y(17);
ro14 = k_dec_Xaa*y(18);
ro15 = k_dec_Xfa*y(19);
ro16 = k_dec_Xc4*y(20);
ro17 = k_dec_Xpro*y(21);
ro18 = k_dec_Xac*y(22)*I18a;
ro18b = k_dec_Xac2*y(37)*I18b; % Acetate Oxidation
ro19 = k_dec_Xh2*y(23);
p_gas_h2 = ( y(33)*R*T_op/16.0 ); % p_gas_h2
p_gas_ch4 =( y(34)*R*T_op/64.0 ); % p_gas_ch4
p_gas_co2 =( y(35)*R*T_op ); % p_gas_co2
p_gas_h2o = (K_H_h2o_base*exp(5290.0*(1.0/T_base - 1.0/T_op)) ); % T
P_gas = ( p_gas_h2+p_gas_ch4+p_gas_co2+p_gas_h2o );
q_gas = k_P*(P_gas-P_atm)*P_gas/P_atm;
if q_gas <0
 q_gas = 0;
end
y1(1) = (q/V_liq)*(S_suf - y(1)) + ro2 +(1-f_fa_li)*ro4 - ro5 ;
% S_aa
y1(2) = (q/V_liq)*(S_aaf - y(2)) + ro3 - ro6;
% S_fa
y1(3) = (q/V_liq)*(S_faf - y(3)) + f_fa_li*ro4 - ro7;
% S_va
y1(4) = ( (q/V_liq)*(S_vaf - y(4)) + (1-Y_aa)*f_va_aa*ro6 - ro8 );
% S_bu
y1(5) = ( (q/V_liq)*(S_buf - y(5)) + (1-Y_su)*f_bu_su*ro5 + (1-Y_aa)*f_bu_aa*ro6 - ro9 );
% S_pro
y1(6) = ( (q/V_liq)*(S_prof - y(6)) + (1-Y_su)*f_pro_su*ro5 + (1-Y_aa)*f_pro_aa*ro6+(1-Y_c4)*0.54*ro8 - ro10 );
% S_ac
y1(7) = ( (q/V_liq)*(S_acf - y(7)) + (1-Y_su)*f_ac_su*ro5 + (1-Y_aa)*f_ac_aa*ro6 + (1-Y_fa)*0.7*ro7 + (1-Y_c4)*0.31*ro8 + (1-Y_c4)*0.8*ro9 + (1-Y_pro)*0.57*ro10 - ro11 -ro11b );
% S_h2
y1(8) = (q/V_liq)*(S_h2f - y(8)) + (1-Y_su)*f_h2_su*ro5 + (1-Y_aa)*f_h2_aa*ro6 + (1-Y_fa)*0.3*ro7 + (1-Y_c4)*0.15*ro8 + (1-Y_c4)*0.2*ro9 +(1-Y_pro)*0.43*ro10 + (1-Y_ac2)*ro11b - ro12 -( kLa*(y(8)-16*K_H_h2*(y(33)*R*T_op/16)) );
% S_ch4
y1(9) = (q/V_liq)*(S_ch4f - y(9)) + (1-Y_ac)*ro11 + (1-Y_h2)*ro12 -( kLa*(y(9)-64*K_H_ch4*(y(34)*R*T_op/64)) );
% S_IC
y1(10) = ( (q/V_liq)*(S_ICf - y(10)) - stoich1*ro1 - stoich2*ro2 - stoich3*ro3 - stoich4*ro4 - stoich5*ro5 - stoich6*ro6 - stoich7*ro7 - stoich8*ro8 - stoich9*ro9 - stoich10*ro10 - stoich11*ro11 - stoich11b*ro11b -stoich12*ro12 - stoich13*ro13-stoich13*ro14-stoich13*ro15-stoich13*ro16-stoich13*ro17- stoich13*ro18-stoich13*ro18b-stoich13*ro19- ( kLa*(y(10)-y(31)- K_H_co2*(y(35)*R*T_op)) ) );
y1(11) = ( (q/V_liq)*(S_INf - y(11))-Y_su*N_bac*ro5+(N_aaY_aa*N_bac)*ro6-Y_fa*N_bac*ro7-Y_c4*N_bac*ro8-Y_c4*N_bac*ro9-Y_pro*N_bac*ro10-Y_ac*N_bac*ro11-Y_ac2*N_bac*ro11bY_h2*N_bac*ro12+(N_bacN_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro18b+ro19)+(N_xc-f_xI_xc*N_If_sI_xc*N_I-f_pr_xc*N_aa)*ro1 );
% S_I
y1(12) = (q/V_liq)*(S_If - y(12))+f_sI_xc*ro1;
% X_c
y1(13) = (q/V_liq)*(X_cf - y(13))- ro1 + ro13 + ro14 + ro15 + ro16 +ro17 + ro18 + ro18b + ro19;
% X_ch
y1(14) = (q/V_liq)*(X_chf - y(14)) + f_ch_xc*ro1 - ro2;
% X_pr
y1(15) = (q/V_liq)*(X_prf - y(15)) + f_pr_xc*ro1 - ro3;
% X_li
y1(16) = (q/V_liq)*(X_lif - y(16))+ f_li_xc*ro1 - ro4;
% X_su
y1(17) = (q/V_liq)*(X_suf - y(17)) + Y_su*ro5 - ro13;
% X_aa
y1(18) = (q/V_liq)*(X_aaf - y(18))+ Y_aa*ro6 - ro14;
% X_fa
y1(19) = (q/V_liq)*(X_faf - y(19)) + Y_fa*ro7 - ro15;
% X_c4
y1(20) = (q/V_liq)*(X_c4f - y(20)) + Y_c4*ro8 + Y_c4*ro9 - ro16;
% X_pro
y1(21) = (q/V_liq)*(X_prof - y(21)) + Y_pro*ro10 - ro17;
% X_ac
y1(22) = (q/V_liq)*(X_acf - y(22)) + Y_ac*ro11 - ro18;
% X_h2
y1(23) = (q/V_liq)*(X_h2f - y(23)) + Y_h2*ro12 - ro19;
% X_I
y1(24) = (q/V_liq)*(X_If - y(24)) + f_xI_xc*ro1;
% S_cat_ion
y1(25) = ( (q/V_liq)*(S_cat_ionf - y(25)) );
% S_an_ion
y1(26) = ( (q/V_liq)*(S_an_ionf - y(26)) );
% S_va_ion
y1(27) = - r_A_4;
% S_bu_ion
y1(28) = - r_A_5;
% S_pro_ion
y1(29) = - r_A_6;
% S_ac_ion
y1(30) = - r_A_7;
% S_hco3_ion
y1(31) = - r_A_10;
% S_nh3_ion
y1(32) = - r_A_11;
% S_gas_h2
y1(33) = - y(33)*(q_gas/V_gas) +r_T_8*(V_liq/V_gas);
% S_gas_ch4
y1(34) = - y(34)*(q_gas/V_gas) + r_T_9*(V_liq/V_gas);
% S_gas_co2
y1(35) = - y(35)*(q_gas/V_gas)+ r_T_10*(V_liq/V_gas);
A1 = ( (q/V_liq)*(S_an_ionf - y(26)) ); % dSan-/dt
A2 = ( (q/V_liq)*(S_INf - y(11))-Y_su*N_bac*ro5+(N_aaY_aa*N_bac)*ro6-Y_fa*N_bac*ro7-Y_c4*N_bac*ro8-Y_c4*N_bac*ro9-Y_pro*N_bac*ro10-Y_ac*N_bac*ro11-Y_h2*N_bac*ro12+(N_bacN_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro19)+(N_xc-f_xI_xc*N_I-f_sI_xc*N_If_pr_xc*N_aa)*ro1 ); % dSIN/dt
A3 = ( (q/V_liq)*(S_ICf - y(10)) - stoich1*ro1-stoich2*ro2-stoich3*ro3-stoich4*ro4-stoich5*ro5-stoich6*ro6-stoich7*ro7-stoich8*ro8-stoich9*ro9-stoich10*ro10-stoich11*ro11-stoich11b*ro11b-stoich12*ro12-stoich13*ro13-stoich13*ro14-stoich13*ro15-stoich13*ro16-stoich13*ro17-stoich13*ro18-stoich13*ro18b-stoich13*ro19- ( kLa*(y(10)-y(31)-K_H_co2*(y(35)*R*T_op)) ) );
A4=( (q/V_liq)*(S_acf - y(7)) + (1-Y_su)*f_ac_su*ro5 + (1-Y_aa)*f_ac_aa*ro6 + (1-Y_fa)*0.7*ro7 + (1-Y_c4)*0.31*ro8 + (1-Y_c4)*0.8*ro9 + (1-Y_pro)*0.57*ro10 - ro11 -ro11b );
A5 = ( (q/V_liq)*(S_prof - y(6)) + (1-Y_su)*f_pro_su*ro5 + (1-Y_aa)*f_pro_aa*ro6+(1-Y_c4)*0.54*ro8 - ro10 );
A6 = ( (q/V_liq)*(S_buf - y(5)) + (1-Y_su)*f_bu_su*ro5 + (1-Y_aa)*f_bu_aa*ro6 - ro9 );
A7 = ( (q/V_liq)*(S_vaf - y(4)) + (1-Y_aa)*f_va_aa*ro6 - ro8 );
A8 = ( (q/V_liq)*(S_INf - y(11))-Y_su*N_bac*ro5+(N_aa-Y_aa*N_bac)*ro6-Y_fa*N_bac*ro7-Y_c4*N_bac*ro8-Y_c4*N_bac*ro9-Y_pro*N_bac*ro10-Y_ac*N_bac*ro11-Y_ac2*N_bac*ro11b-Y_h2*N_bac*ro12+(N_bacN_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro18b+ro19)+(N_xc-f_xI_xc*N_If_sI_xc*N_I-f_pr_xc*N_aa)*ro1 );
A9= ( (q/V_liq)*(S_cat_ionf - y(25)) );
A =A1+A2*K_a_IN/(K_a_IN+y(36))+A3*K_a_co2/(K_a_co2+y(36))+(1/64)*A4*K_a_ac/(K_a_ac+y(36)) + (1/112)*A5*K_a_pro/(K_a_pro+y(36)) +(1/160)*A6*K_a_bu/(K_a_bu+y(36)) + (1/208)*A7*K_a_va/(K_a_va+y(36)) -A8 - A9;
B = 1 + y(11)*K_a_IN/((K_a_IN+y(36))^2) +y(10)*K_a_co2/((K_a_co2+y(36))^2) +(1/64)*y(7)*K_a_ac/((K_a_ac+y(36))^2) +(1/112)*y(6)*K_a_pro/((K_a_pro+y(36))^2) +(1/160)*y(5)*K_a_bu/((K_a_bu+y(36))^2) +(1/208)*y(4)*K_a_va/((K_a_va+y(36))^2) + K_w/(y(36)^2);
y1(36) = A/B; 
y1(37) = (q/V_liq)*(X_ac2f - y(37)) + Y_ac2*ro11b - ro18b;
clear q_gas

And That's for the ODE solving and plots (AD.m)

function ad(tspan,u)
format long
[t,y] = ode15s(ADDE,tspan,u);
global Y_ac2;
global k_m_ac2
global k_m_su;
global K_S_su;
global k_m_aa;
global K_S_aa;
global k_m_fa;
global K_S_fa;
global k_m_c4;
global K_S_c4;
global k_m_pro;
global K_S_pro;
global k_m_ac;
global K_S_ac;
global k_m_h2;
global K_S_h2;
global inhib11b;
global inhib11;
global inhib12;
global inhib56;
global inhib7;
global inhib10;
global Y_su;
global f_h2_su;
global Y_aa;
global f_h2_aa;
global Y_fa;
global Y_c4;
global Y_pro;
global Y_ac;
global Y_h2;
global minx
global R
global T_op
global K_H_ch4
global K_H_co2
global K_H_h2
global kLa
global maxx
global K_S_ac2
global af
af = 1.2; % adjustment factor for accurate concentration prediction
minx = 0;
figure (1)
set(gcf, 'color', [1 1 1])
subplot(3,2,1);
plot(t,(T_op/273.15)*1.4*( kLa*(y(:,8)-16*K_H_h2*(y(:,33)*R*T_op/16)) ),'Linewidth',2,'color','r','LineStyle','-');
title('Volumetric production of gas H_2','FontSize',10,'Fontname','cmr10');
xlabel(['Time ~ ',num2str(tspan(1,2),'%5.4g'),' (days)'],'FontSize',10,'Fontname','cmr10'), ylabel('m^3 H_2 m^-^3 d^-^1','Rotation',90,'FontSize',10,'Fontname','cmr10'),
xlim([minx maxx])
set(gca,'fontsize',10,'Fontname','cmr10');
subplot(3,2,2);
plot(t,(T_op/273.15)*0.35*( kLa*(y(:,9)-64*K_H_ch4*(y(:,34)*R*T_op/64)) ),'Linewidth',2,'color','g','LineStyle','-');
title('Volumetric production of gas CH_4','FontSize',10,'Fontname','cmr10');
xlabel(['Time ~ ',num2str(tspan(1,2),'%5.4g'),' (days)'],'FontSize',10,'Fontname','cmr10'), ylabel('m^3 CH_4 m^-^3 d^-^1','Rotation',90,'FontSize',10, 'Fontname','cmr10'),
xlim([minx maxx])
set(gca,'fontsize',10,'Fontname','cmr10');
subplot(3,2,3);
plot(t,(T_op/273.15)*22.4*( kLa*(y(:,10)-y(:,31)-K_H_co2*(y(:,35)*R*T_op)) ) ,'Linewidth',2, 'color','b','LineStyle','-');
title('Volumetric production of gas CO_2','FontSize',10,'Fontname','cmr10');
xlabel(['Time ~ ',num2str(tspan(1,2),'%5.4g'),' (days)'],'FontSize',10,'Fontname','cmr10'), ylabel('m^3 CO_2 m^-^3 d^-^1','Rotation',90,'FontSize',10,'Fontname','cmr10'),xlim([minx maxx])
set(gca,'fontsize',10,'Fontname','cmr10');
subplot(3,2,4);
plot(t,(T_op/273.15)*1.4*( kLa*(y(:,8)-16*K_H_h2*(y(:,33)*R*T_op/16)) )+(T_op/273.15)*0.35*( kLa*(y(:,9)-64*K_H_ch4*(y(:,34)*R*T_op/64)) )+(T_op/273.15)*22.4*( kLa*(y(:,10)-y(:,31)-K_H_co2*(y(:,35)*R*T_op)) ) ,'Linewidth',2,'color','b','LineStyle','-');
xlim([minx maxx])
title('Volumetric production of Biogas','FontSize',10,'Fontname','cmr10');
xlabel(['Time ~ ',num2str(tspan(1,2),'%5.4g'),' (days)'],...
 'FontSize',10,'Fontname','cmr10'), ylabel('m^3 biogas m^-^3 d^-^1','Rotation',90,'FontSize',10,'Fontname','cmr10'),
set(gca,'fontsize',10,'Fontname','cmr10');
subplot(3,1,3);
plot(t, ( T_op/273.15)*0.35*( kLa*(y(:,9)-64*K_H_ch4*(y(:,34)*R*T_op/64)) ) ./( (T_op/273.15)*1.4*( kLa*(y(:,8)-16*K_H_h2*(y(:,33)*R*T_op/16)) )+(T_op/273.15)*0.35*( kLa*(y(:,9)-64*K_H_ch4*(y(:,34)*R*T_op/64)) )+(T_op/273.15)*22.4*( kLa*(y(:,10)-y(:,31)-K_H_co2*(y(:,35)*R*T_op)) ) ) ,'Linewidth',2,'color','g','LineStyle','-');
hold on
plot(t, (T_op/273.15)*22.4*( kLa*(y(:,10)-y(:,31)-K_H_co2*(y(:,35)*R*T_op)) ) ./ ( (T_op/273.15)*1.4*( kLa*(y(:,8)-16*K_H_h2*(y(:,33)*R*T_op/16)) )+(T_op/273.15)*0.35*( kLa*(y(:,9)-64*K_H_ch4*(y(:,34)*R*T_op/64)) )+(T_op/273.15)*22.4*( kLa*(y(:,10)-y(:,31)-K_H_co2*(y(:,35)*R*T_op)) ) ) ,'Linewidth',2,'color','b','LineStyle','-');
xlim([minx maxx])
title('% of CH_4 and CO_2 in biogas','FontSize',10, 'Fontname','cmr10');
xlabel(['Time ~ ',num2str(tspan(1,2),'%5.4g'),' (days)'],...
 'FontSize',10,'Fontname','cmr10'), ylabel('% inBiogas)','Rotation',90,'FontSize',10,'Fontname','cmr10'),
hold off
set(gca,'fontsize',10,'Fontname','cmr10');
legend('CH_4','CO_2')
figure(2)
set(gcf, 'color', [1 1 1])
subplot(3,2,1);
plot(t,af*(1000*y(:,4)),'Linewidth',2, 'color','r','LineStyle','-');
title('Valeric and Butyric Acids','FontSize',10, 'Fontname','cmr10');
xlabel(['Time = ',num2str(maxx,'%5.4g'),' (days)'],'FontSize',10,'Fontname','cmr10'), ylabel('Valeric & Butyric mg L^-^1','Rotation',90,'FontSize',10,'Fontname','cmr10'),
xlim([minx maxx])
hold on
plot(t,af*(1000*y(:,5)),'Linewidth',2, 'color','c','LineStyle','-.');
xlim([minx maxx])
hold off
legend('Valeric','Butyric')
set(gca,'fontsize',10,'Fontname','cmr10');
subplot(3,2,2);
plot(t,af*(1000*y(:,6)),'Linewidth',2, 'color','b','LineStyle','-');
title('Propionic Acid','FontSize',10, 'Fontname','cmr10');
xlabel(['Time =',num2str(maxx,'%5.4g'),' (days)'],'FontSize',10,'Fontname','cmr10'), ylabel('Propionic mg L^-^1','Rotation',90,'FontSize',10,'Fontname','cmr10'),
xlim([minx maxx])
set(gca,'fontsize',10,'Fontname','cmr10');
subplot(3,2,3);
plot(t,af*(1000*y(:,7)),'Linewidth',2, 'color','m','LineStyle','-');
title('Acetic Acid','FontSize',10, 'Fontname','cmr10');
xlabel(['Time = ',num2str(maxx,'%5.4g'),' (days)'],'FontSize',10,'Fontname','cmr10'), ylabel('Acetic mg L^-^1','Rotation',90,'FontSize',10,'Fontname','cmr10'),
xlim([minx maxx])
set(gca,'fontsize',10,'Fontname','cmr10');
subplot(3,2,4);
% convert kg COD/ m3 to mg /L
plot(t,af*(1000*(y(:,4)+y(:,5)+y(:,6)+y(:,7))),'Linewidth',2,'color','k','LineStyle','-');
title('Total Volatile Fatty Acid','FontSize',10, 'Fontname','cmr10');
xlabel(['Time = ',num2str(maxx,'%5.4g'),' (days)'],'FontSize',10,'Fontname','cmr10'), ylabel('Total VFAs (mg L^-^1)','Rotation',90,'FontSize',10,'Fontname','cmr10'),
xlim([minx maxx])
set(gca,'fontsize',10,'Fontname','cmr10');

Many thanks for your time and attached is the excel file for the inputs

[Walter Roberson]

First you need to go to

function ad(tspan,u)
format long
[t,y] = ode15s(ADDE,tspan,u);

and change that to

function ad(tspan,u)
format long
[t,y] = ode15s(@ADDE,tspan,u);

after that, run the function that sets the global variables.

Once the global variables have been set, at the command line invoke

global tspan u
ad(tspan,u)

[Mohamed Eisa]

Thanks a lot. In the ADDE code that considers all the calculations, I am finding a new error with it now (not enough inputs) but in line 182, when it starts to take into account y values as y is not defined before in the code. what do you think about it; here is the code:

function y1 = ADDE(t,y)

format long
global q
global V_liq % Volume of liquid part
global V_gas % Volume of gas space
global S_suf
global S_aaf
global S_faf
global S_vaf
global S_buf
global S_prof
global S_acf
global S_h2f
global S_ch4f
global S_ICf
global S_INf
global S_If
global X_cf
global X_chf
global X_prf
global X_lif
global X_suf
global X_aaf
global X_faf
global X_c4f
global X_prof
global X_acf
global X_ac2f
global X_h2f
global X_If
global S_cat_ionf
global S_an_ionf
global f_sI_xc
global f_xI_xc
global f_ch_xc
global f_pr_xc
global f_li_xc
global f_fa_li
global f_h2_su
global f_bu_su
global f_pro_su
global f_ac_su
global f_h2_aa
global f_va_aa
global f_bu_aa
global f_pro_aa
global f_ac_aa
global C_xc
global C_sI
global C_ch
global C_pr
global C_li
global C_xI
global C_su
global C_aa
global C_fa
global C_bu
global C_pro
global C_ac
global C_bac
global C_va
global C_ch4
global N_xc
global N_I
global N_aa
global N_bac
global Y_su
global Y_aa
global Y_fa
global Y_c4
global Y_pro
global Y_ac
global Y_ac2
global Y_h2
global k_dis
global k_hyd_ch
global k_hyd_pr
global k_hyd_li
global k_m_su
global k_m_aa
global k_m_fa
global k_m_c4
global k_m_pro
global k_m_ac
global k_m_ac2
global k_m_h2
global k_dec_Xsu
global k_dec_Xaa
global k_dec_Xfa
global k_dec_Xc4
global k_dec_Xpro
global k_dec_Xac
global k_dec_Xac2
global k_dec_Xh2
global K_S_IN
global K_S_su
global K_S_aa
global K_S_fa
global K_Ih2_fa
global K_S_pro
global K_Ih2_pro
global K_S_ac
global K_S_ac2
global K_I_nh3
global K_S_c4
global K_S_h2
global K_Ih2_c4
global K_Ih2_ac
global kLa
global K_H_h2o_base
global K_H_co2_base
global K_H_ch4_base
global K_H_h2_base
global k_P
global P_atm
global T_base
global T_op
global R
global pK_w_base
global pK_a_va_base
global pK_a_bu_base
global pK_a_pro_base
global pK_a_ac_base
global pK_a_co2_base
global pK_a_IN_base
global k_A_Bva
global k_A_Bbu
global k_A_Bpro
global k_A_Bac
global k_A_Bco2
global k_A_BIN
global pHLim_aa;
global pHLim_ac;
global pHLim_ac2;
global pHLim_h2;
global k_aa;
global k_ac;
global k_ac2;
global k_h2;
global I11a;
global I11b;
global I18a;
global I18b;
global inhib11b;
global inhib56;
global inhib7;
global inhib89;
global inhib10;
global inhib11;
global inhib12;
global Itec4;
global Itepro;
global K_a_va
global K_a_bu
global K_a_pro
global K_a_ac
global K_a_co2
global K_a_IN
global K_w
global K_H_h2
global K_H_ch4
global K_H_co2
global p_gas_h2o
global factor
global P_gas
global p_gas_h2
global p_gas_ch4
global p_gas_co2
global q_gas
factor = (1.0/T_base - 1.0/T_op)/(100*R);
K_a_va =10^(-pK_a_va_base);
K_a_bu = 10^(-pK_a_bu_base);
K_a_pro = 10^(-pK_a_pro_base);
K_a_ac = 10^(-pK_a_ac_base);
K_a_co2 = (10^(-pK_a_co2_base))*exp(7646.0*factor);
K_a_IN = 10^(-pK_a_IN_base)*exp(51965.0*factor); 
K_w = 10^(-pK_w_base)*exp(55900.0*factor); % T adjustment for K_w
K_H_h2 = K_H_h2_base*exp(-4180.0*factor); %/* T adjustment for K_H_h2
K_H_ch4 = K_H_ch4_base*exp(-14240.0*factor);
K_H_co2 = K_H_co2_base*exp(-19410.0*factor);
r_A_4 = ( k_A_Bva*(y(27)*(K_a_va + y(36)) - K_a_va*y(4)) );
r_A_5 = ( k_A_Bbu*(y(28)*(K_a_bu + y(36)) - K_a_bu*y(5)) );
r_A_6 = ( k_A_Bpro*(y(29)*(K_a_pro + y(36))- K_a_pro*y(6)) );
r_A_7 = ( k_A_Bac*(y(30)*(K_a_ac + y(36)) - K_a_ac*y(7)) );
r_A_10 =( k_A_Bco2*(y(31)*(K_a_co2 + y(36)) - K_a_co2*y(10)) ); %
r_A_11 =( k_A_BIN*(y(32)*(K_a_IN + y(36)) - K_a_IN*y(11)) ); %
r_T_8 = ( kLa*(y(8)-16*K_H_h2*( y(33)*R*T_op/16.0 )) );
r_T_9 = ( kLa*(y(9)-64*K_H_ch4*( y(34)*R*T_op/64.0 )) );
r_T_10 =( kLa*(y(10)-y(31)-K_H_co2*( y(35)*R*T_op )) );
stoich1 = -C_xc+f_sI_xc*C_sI+f_ch_xc*C_ch+f_pr_xc*C_pr+f_li_xc*C_li+f_xI_xc*C_xI;
stoich2 = -C_ch+C_su;
stoich3 = -C_pr+C_aa;
stoich4 = -C_li+(1.0-f_fa_li)*C_su+f_fa_li*C_fa;
stoich5 = -C_su+(1.0-Y_su)*(f_bu_su*C_bu+f_pro_su*C_pro+f_ac_su*C_ac)+Y_su*C_bac;
stoich6 = -C_aa+(1.0-Y_aa)*(f_va_aa*C_va+f_bu_aa*C_bu+f_pro_aa*C_pro+f_ac_aa*C_ac)+Y_aa*C_bac;
stoich7 = -C_fa+(1.0-Y_fa)*0.7*C_ac+Y_fa*C_bac;
stoich8 = -C_va+(1.0-Y_c4)*0.54*C_pro+(1.0-Y_c4)*0.31*C_ac+Y_c4*C_bac;
stoich9 = -C_bu+(1.0-Y_c4)*0.8*C_ac+Y_c4*C_bac;
stoich10 = -C_pro+(1.0-Y_pro)*0.57*C_ac+Y_pro*C_bac;
stoich11 = -C_ac+(1.0-Y_ac)*C_ch4+Y_ac*C_bac;
stoich11b = -C_ac+Y_ac2*C_bac;
stoich12 = (1.0-Y_h2)*C_ch4+Y_h2*C_bac;
stoich13 = -C_bac+C_xc;
I_pH_aa = pHLim_aa^k_aa /(y(36)^k_aa + pHLim_aa^k_aa);
I_pH_ac = pHLim_ac^k_ac /(y(36)^k_ac + pHLim_ac^k_ac);
I_pH_ac2 = pHLim_ac2^k_ac2 /(y(36)^k_ac2 + pHLim_ac2^k_ac2);
I_pH_h2 = pHLim_h2^k_h2 /(y(36)^k_h2 + pHLim_h2^k_h2);
I_IN_lim = 1.0/(1.0+K_S_IN/y(11));
I_h2_fa = 1.0/(1.0+y(8)/K_Ih2_fa);
I_h2_c4 = 1.0/(1.0+y(8)/K_Ih2_c4);
I_h2_pro = 1.0/(1.0+y(8)/K_Ih2_pro);
I_h2_ac = 1.0/(1.0+y(8)/K_Ih2_ac);
I_nh3 = 1.0/(1.0+y(32)/K_I_nh3);
Itec4 = 1;
Itepro = 1;
inhib56 = I_pH_aa*I_IN_lim;
inhib7 = inhib56*I_h2_fa;
inhib89 = inhib56*I_h2_c4*Itec4;
inhib10 = inhib56*I_h2_pro*Itepro;
inhib11 = I_pH_ac*I_IN_lim*I_nh3*I11a;
inhib11b = I_pH_ac2*I_IN_lim*I_h2_ac*I11b;
inhib12 = I_pH_h2*I_IN_lim;
ro1 = k_dis*y(13);
ro2 = k_hyd_ch*y(14);
ro3 = k_hyd_pr*y(15);
ro4 = k_hyd_li*y(16);
ro5 = k_m_su*(y(1)/(y(1)+K_S_su))*y(17)*inhib56;
ro6 = k_m_aa*(y(2)/(K_S_aa+y(2)))*y(18)*inhib56;
ro7 = k_m_fa*(y(3)/(K_S_fa+y(3)))*y(19)*inhib7;
ro8 = k_m_c4*(y(4)/(K_S_c4+y(4)))*y(20)*(y(4)/(y(5)+y(4)+1e-6))*inhib89;
ro9 = k_m_c4*(y(5)/(K_S_c4+y(5)))*y(20)*(y(5)/(y(4)+y(5)+1e-6))*inhib89;
ro10 = k_m_pro*(y(6)/(K_S_pro+y(6)))*y(21)*inhib10;
ro11 = k_m_ac*(y(7)/(K_S_ac+y(7)))*y(22)*inhib11;
ro11b = k_m_ac2*(y(7)/(K_S_ac2+y(7)))*y(37)*inhib11b; % Acetate Oxidation
ro12 = k_m_h2*(y(8)/(K_S_h2+y(8)))*y(23)*inhib12;
ro13 = k_dec_Xsu*y(17);
ro14 = k_dec_Xaa*y(18);
ro15 = k_dec_Xfa*y(19);
ro16 = k_dec_Xc4*y(20);
ro17 = k_dec_Xpro*y(21);
ro18 = k_dec_Xac*y(22)*I18a;
ro18b = k_dec_Xac2*y(37)*I18b; % Acetate Oxidation
ro19 = k_dec_Xh2*y(23);
p_gas_h2 = ( y(33)*R*T_op/16.0 ); % p_gas_h2
p_gas_ch4 =( y(34)*R*T_op/64.0 ); % p_gas_ch4
p_gas_co2 =( y(35)*R*T_op ); % p_gas_co2
p_gas_h2o = (K_H_h2o_base*exp(5290.0*(1.0/T_base - 1.0/T_op)) ); % T
P_gas = ( p_gas_h2+p_gas_ch4+p_gas_co2+p_gas_h2o );
q_gas = k_P*(P_gas-P_atm)*P_gas/P_atm;
if q_gas <0
 q_gas = 0;
end
y1(1) = (q/V_liq)*(S_suf - y(1)) + ro2 +(1-f_fa_li)*ro4 - ro5 ;
% S_aa
y1(2) = (q/V_liq)*(S_aaf - y(2)) + ro3 - ro6;
% S_fa
y1(3) = (q/V_liq)*(S_faf - y(3)) + f_fa_li*ro4 - ro7;
% S_va
y1(4) = ( (q/V_liq)*(S_vaf - y(4)) + (1-Y_aa)*f_va_aa*ro6 - ro8 );
% S_bu
y1(5) = ( (q/V_liq)*(S_buf - y(5)) + (1-Y_su)*f_bu_su*ro5 + (1-Y_aa)*f_bu_aa*ro6 - ro9 );
% S_pro
y1(6) = ( (q/V_liq)*(S_prof - y(6)) + (1-Y_su)*f_pro_su*ro5 + (1-Y_aa)*f_pro_aa*ro6+(1-Y_c4)*0.54*ro8 - ro10 );
% S_ac
y1(7) = ( (q/V_liq)*(S_acf - y(7)) + (1-Y_su)*f_ac_su*ro5 + (1-Y_aa)*f_ac_aa*ro6 + (1-Y_fa)*0.7*ro7 + (1-Y_c4)*0.31*ro8 + (1-Y_c4)*0.8*ro9 + (1-Y_pro)*0.57*ro10 - ro11 -ro11b );
% S_h2
y1(8) = (q/V_liq)*(S_h2f - y(8)) + (1-Y_su)*f_h2_su*ro5 + (1-Y_aa)*f_h2_aa*ro6 + (1-Y_fa)*0.3*ro7 + (1-Y_c4)*0.15*ro8 + (1-Y_c4)*0.2*ro9 +(1-Y_pro)*0.43*ro10 + (1-Y_ac2)*ro11b - ro12 -( kLa*(y(8)-16*K_H_h2*(y(33)*R*T_op/16)) );
% S_ch4
y1(9) = (q/V_liq)*(S_ch4f - y(9)) + (1-Y_ac)*ro11 + (1-Y_h2)*ro12 -( kLa*(y(9)-64*K_H_ch4*(y(34)*R*T_op/64)) );
% S_IC
y1(10) = ( (q/V_liq)*(S_ICf - y(10)) - stoich1*ro1 - stoich2*ro2 - stoich3*ro3 - stoich4*ro4 - stoich5*ro5 - stoich6*ro6 - stoich7*ro7 - stoich8*ro8 - stoich9*ro9 - stoich10*ro10 - stoich11*ro11 - stoich11b*ro11b -stoich12*ro12 - stoich13*ro13-stoich13*ro14-stoich13*ro15-stoich13*ro16-stoich13*ro17- stoich13*ro18-stoich13*ro18b-stoich13*ro19- ( kLa*(y(10)-y(31)- K_H_co2*(y(35)*R*T_op)) ) );
y1(11) = ( (q/V_liq)*(S_INf - y(11))-Y_su*N_bac*ro5+(N_aaY_aa*N_bac)*ro6-Y_fa*N_bac*ro7-Y_c4*N_bac*ro8-Y_c4*N_bac*ro9-Y_pro*N_bac*ro10-Y_ac*N_bac*ro11-Y_ac2*N_bac*ro11bY_h2*N_bac*ro12+(N_bacN_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro18b+ro19)+(N_xc-f_xI_xc*N_If_sI_xc*N_I-f_pr_xc*N_aa)*ro1 );
% S_I
y1(12) = (q/V_liq)*(S_If - y(12))+f_sI_xc*ro1;
% X_c
y1(13) = (q/V_liq)*(X_cf - y(13))- ro1 + ro13 + ro14 + ro15 + ro16 +ro17 + ro18 + ro18b + ro19;
% X_ch
y1(14) = (q/V_liq)*(X_chf - y(14)) + f_ch_xc*ro1 - ro2;
% X_pr
y1(15) = (q/V_liq)*(X_prf - y(15)) + f_pr_xc*ro1 - ro3;
% X_li
y1(16) = (q/V_liq)*(X_lif - y(16))+ f_li_xc*ro1 - ro4;
% X_su
y1(17) = (q/V_liq)*(X_suf - y(17)) + Y_su*ro5 - ro13;
% X_aa
y1(18) = (q/V_liq)*(X_aaf - y(18))+ Y_aa*ro6 - ro14;
% X_fa
y1(19) = (q/V_liq)*(X_faf - y(19)) + Y_fa*ro7 - ro15;
% X_c4
y1(20) = (q/V_liq)*(X_c4f - y(20)) + Y_c4*ro8 + Y_c4*ro9 - ro16;
% X_pro
y1(21) = (q/V_liq)*(X_prof - y(21)) + Y_pro*ro10 - ro17;
% X_ac
y1(22) = (q/V_liq)*(X_acf - y(22)) + Y_ac*ro11 - ro18;
% X_h2
y1(23) = (q/V_liq)*(X_h2f - y(23)) + Y_h2*ro12 - ro19;
% X_I
y1(24) = (q/V_liq)*(X_If - y(24)) + f_xI_xc*ro1;
% S_cat_ion
y1(25) = ( (q/V_liq)*(S_cat_ionf - y(25)) );
% S_an_ion
y1(26) = ( (q/V_liq)*(S_an_ionf - y(26)) );
% S_va_ion
y1(27) = - r_A_4;
% S_bu_ion
y1(28) = - r_A_5;
% S_pro_ion
y1(29) = - r_A_6;
% S_ac_ion
y1(30) = - r_A_7;
% S_hco3_ion
y1(31) = - r_A_10;
% S_nh3_ion
y1(32) = - r_A_11;
% S_gas_h2
y1(33) = - y(33)*(q_gas/V_gas) +r_T_8*(V_liq/V_gas);
% S_gas_ch4
y1(34) = - y(34)*(q_gas/V_gas) + r_T_9*(V_liq/V_gas);
% S_gas_co2
y1(35) = - y(35)*(q_gas/V_gas)+ r_T_10*(V_liq/V_gas);
A1 = ( (q/V_liq)*(S_an_ionf - y(26)) ); % dSan-/dt
A2 = ( (q/V_liq)*(S_INf - y(11))-Y_su*N_bac*ro5+(N_aaY_aa*N_bac)*ro6-Y_fa*N_bac*ro7-Y_c4*N_bac*ro8-Y_c4*N_bac*ro9-Y_pro*N_bac*ro10-Y_ac*N_bac*ro11-Y_h2*N_bac*ro12+(N_bacN_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro19)+(N_xc-f_xI_xc*N_I-f_sI_xc*N_If_pr_xc*N_aa)*ro1 ); % dSIN/dt
A3 = ( (q/V_liq)*(S_ICf - y(10)) - stoich1*ro1-stoich2*ro2-stoich3*ro3-stoich4*ro4-stoich5*ro5-stoich6*ro6-stoich7*ro7-stoich8*ro8-stoich9*ro9-stoich10*ro10-stoich11*ro11-stoich11b*ro11b-stoich12*ro12-stoich13*ro13-stoich13*ro14-stoich13*ro15-stoich13*ro16-stoich13*ro17-stoich13*ro18-stoich13*ro18b-stoich13*ro19- ( kLa*(y(10)-y(31)-K_H_co2*(y(35)*R*T_op)) ) );
A4=( (q/V_liq)*(S_acf - y(7)) + (1-Y_su)*f_ac_su*ro5 + (1-Y_aa)*f_ac_aa*ro6 + (1-Y_fa)*0.7*ro7 + (1-Y_c4)*0.31*ro8 + (1-Y_c4)*0.8*ro9 + (1-Y_pro)*0.57*ro10 - ro11 -ro11b );
A5 = ( (q/V_liq)*(S_prof - y(6)) + (1-Y_su)*f_pro_su*ro5 + (1-Y_aa)*f_pro_aa*ro6+(1-Y_c4)*0.54*ro8 - ro10 );
A6 = ( (q/V_liq)*(S_buf - y(5)) + (1-Y_su)*f_bu_su*ro5 + (1-Y_aa)*f_bu_aa*ro6 - ro9 );
A7 = ( (q/V_liq)*(S_vaf - y(4)) + (1-Y_aa)*f_va_aa*ro6 - ro8 );
A8 = ( (q/V_liq)*(S_INf - y(11))-Y_su*N_bac*ro5+(N_aa-Y_aa*N_bac)*ro6-Y_fa*N_bac*ro7-Y_c4*N_bac*ro8-Y_c4*N_bac*ro9-Y_pro*N_bac*ro10-Y_ac*N_bac*ro11-Y_ac2*N_bac*ro11b-Y_h2*N_bac*ro12+(N_bacN_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro18b+ro19)+(N_xc-f_xI_xc*N_If_sI_xc*N_I-f_pr_xc*N_aa)*ro1 );
A9= ( (q/V_liq)*(S_cat_ionf - y(25)) );
A =A1+A2*K_a_IN/(K_a_IN+y(36))+A3*K_a_co2/(K_a_co2+y(36))+(1/64)*A4*K_a_ac/(K_a_ac+y(36)) + (1/112)*A5*K_a_pro/(K_a_pro+y(36)) +(1/160)*A6*K_a_bu/(K_a_bu+y(36)) + (1/208)*A7*K_a_va/(K_a_va+y(36)) -A8 - A9;
B = 1 + y(11)*K_a_IN/((K_a_IN+y(36))^2) +y(10)*K_a_co2/((K_a_co2+y(36))^2) +(1/64)*y(7)*K_a_ac/((K_a_ac+y(36))^2) +(1/112)*y(6)*K_a_pro/((K_a_pro+y(36))^2) +(1/160)*y(5)*K_a_bu/((K_a_bu+y(36))^2) +(1/208)*y(4)*K_a_va/((K_a_va+y(36))^2) + K_w/(y(36)^2);
y1(36) = A/B; 
y1(37) = (q/V_liq)*(X_ac2f - y(37)) + Y_ac2*ro11b - ro18b;
clear q_gas

[Walter Roberson]

There are two references to the undefined variable N_aaY_aa

(There are other problems that I spent some time cleaning up, but this one I cannot deduce the value of.)

[Mohamed Eisa]

Thanks for your efforts, sir! I can't thank you enough! if it is possible if you could share with me the edited code by you, that would be totally appreciated.

Yours,

Mohamed

[Walter Roberson]

I would prefer if you posted the definition of N_aaY_aa so I can do more clean up and testing.