Model ADM1, using ODEs

I'm trying to model ADM1 in matlab, however my graphics generated when running show errors during execution, I would like to understand what is happening, my codes are divided into 3 parts: globalvariables, ADDE, ad.

The second has the derivatives that will be solved and in the third it calls ADDE to solve.

globalvariables

format long

clearvars
%---------------------------
%---------------------------
% Digester configuration and tspan
global q V_dig V_liq V_gas tspan u maxx

%---------------------------
%---------------------------
% Variables of soluble and particulate components come in digester
global S_suf S_aaf S_faf S_vaf S_buf S_prof S_acf S_h2f S_ch4f S_ICf S_INf S_If 
global X_cf X_chf X_prf X_lif X_suf X_aaf X_faf X_c4f X_prof X_acf X_ac2f X_h2f X_If
global S_cat_ionf S_an_ionf S_va_ionf S_bu_ionf S_pro_ionf S_ac_ionf S_hco3_ionf S_nh3f S_gas_h2f S_gas_ch4f S_gas_co2f S_h_ionf

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%---------------------------------
%Fraction of each components in Xc
global f_sI_xc f_xI_xc f_ch_xc f_pr_xc f_li_xc f_fa_li f_h2_su f_bu_su f_pro_su f_ac_su f_h2_aa f_va_aa f_bu_aa f_pro_aa f_ac_aa

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% -----------------------------------------------
% Carbon and Nitrogen concentration in components
global C_xc C_sI C_ch C_pr C_li C_xI C_su C_aa C_fa C_bu C_pro C_ac C_bac C_va C_ch4
global N_xc N_I N_aa N_bac

%^^^^^^^^^^^^^^^^^^^^^^^^
%------------------------
% Yield uptake Components
global Y_su Y_aa Y_fa Y_c4 Y_pro Y_ac Y_ac2 Y_h2

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%-----------------------------------------------------
% RATES OF disintegration, hydrolysis and coefficients
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

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%-----------------------------
% Half saturation coefficients
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

%^^^^^^^^^^^^^^^^^^^^^^^^
%------------------------
% Acid and Gas parameters
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


%^^^^^^^^^^^^^^^^^^^^^^^^
%------------------------
% Acid and Gas parameters
global kLa K_H_h2o_base K_H_co2_base K_H_ch4_bas K_H_h2_base k_P
global P_atm
global T_base T_op
global R
global pK_w_base pK_a_va_base pK_a_bu_base pK_a_pro_base pK_a_ac_base pK_a_co2_base pK_a_IN_base k_A_Bva k_A_Bbu k_A_Bpro k_A_Bac k_A_Bco2 k_A_BIN
global pH_UL_h2 pH_LL_h2 pH_UL_aa pH_LL_aa pH_UL_ac pH_LL_ac pH_UL_ac2 pH_LL_ac2
global K_Ih2_ac

%^^^^^^^^^^^^^^^^^^^^^^^^
%------------------------
% Inhibition factors
global pHLim_aa pHLim_ac pHLim_ac2 pHLim_h2
global k_aa k_ac k_ac2 k_h2
global I11a I11b I18a I18b


%
% ------------------------------------------------------------------------
%-------- Read input data from excel file
u = xlsread('ADM1data.xlsx','Digesterconfig','K3:K39'); % Initial conditions for DEs
Inputs = xlsread('ADM1data.xlsx','Digesterconfig','F3:F39');
Digesterconfig = xlsread('ADM1data.xlsx','Digesterconfig','B3:B9');
Fraction = xlsread('ADM1data.xlsx','Biostoic','C3:C17');
Carbonstoichiometries = xlsread('ADM1data.xlsx','Biostoic','C21:C35');
Nitrogenstoichiometries = xlsread('ADM1data.xlsx','Biostoic','C39:C42');
Yielduptakecomponents = xlsread('ADM1data.xlsx','Biochemrate','C3:C10');
Dishydcoefficients = xlsread('ADM1data.xlsx','Biochemrate','C13:C32');
Halfsaturatecoefficients = xlsread('ADM1data.xlsx','Biochemrate','C34:C47');
Acidgasparameters = xlsread('ADM1data.xlsx','Biochemrate','C50:C80');

%
% Assign values to variables
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);

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%----------------------------------------------------------------------
% The method suggested by Siegrist et al. (2002) used a Hill inhibition
% function based on the hydrogen ion concentration instead to calculate
% inhibition factors.
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 = 24/(pH_UL_aa-pH_LL_aa);
k_ac = 45/(pH_UL_ac-pH_LL_ac);
k_ac2 = 45/(pH_UL_ac2-pH_LL_ac2);
k_h2 = 3/(pH_UL_h2-pH_LL_h2);

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%----------------------------------------------------------
% Setup initial condition for running both pathways AC & AO
%Para rodar o código para o ADM1 original I11b e I18b deve ser =0
I11a = 1;
I11b = 0;
I18a = 1;
I18b = 0;

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% run the model (alterar depois que o código esteja finalizado
%ad(tspan,u)

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%-----------------------------------
% Delete all the temporary variables
clear Inputs;
clear Digesterconfig;
clear Fraction;
clear Carbonstoichiometries;
clear Nitrogenstoichiometries;
clear Yielduptakecomponents;
clear Dishydcoefficients;
clear Halfsaturatecoefficients;
clear Acidgasparameters;

%
% sabe the results
%save saveddata
save('results.mat')

ADDE

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

format long

% Encontra-se as derivadas
% O método utilizado para resolver as dt é o método de Euler
% Adiciono o produto de um tamanho e ocorrem mudanças nas variaveis 

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%----------------------------------
% Digester configurations and tspan

global q % Flow
global V_liq % Volume of liquid part
global V_gas % Volume of gas space

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%-----------------------------------------------------------------
% Variables of soluble and particulate components come in digester

global S_suf S_aaf S_faf S_vaf S_buf S_prof S_acf S_h2f S_ch4f S_ICf S_INf S_If
global X_cf X_chf X_prf X_lif X_suf X_aaf X_faf X_c4f X_prof X_acf X_ac2f X_h2f X_If
global S_cat_ionf S_an_ionf

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%---------------------------------
%Fraction of each components in Xc
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

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% -----------------------------------------------
% Carbon and Nitrogen concentration in components
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

%^^^^^^^^^^^^^^^^^^^^^^^^
%------------------------
% Yield uptake Components
global Y_su
global Y_aa
global Y_fa
global Y_c4
global Y_pro
global Y_ac
global Y_ac2
global Y_h2

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%-----------------------------------------------------
% RATES OF disintegration, hydrolysis and coefficients
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

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%-----------------------------
% Half saturation coefficients
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

%^^^^^^^^^^^^^^^^^^^^^^^^
%------------------------
% Acid and Gas parameters
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

%^^^^^^^^^^^^^^^^^^^^^^^^
%------------------------
% Inhibition factors
global I11a % - pHac + INlim + H2ac
global I11b % - pHac + INlim + H2ac2
global I18a %decay of Xac
global I18b %decay of Xac
global inhib11b
global inhib56 %su and aa
global inhib7 %LCFA
global inhib89 %va and bu
global inhib10 %pro
global inhib11 %ac
global inhib12 %h2
global Itec4 %inibição de fatores traços de va e bu uptake
global Itepro %inibição de fatores traços de pro uptake

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%--------------------------------------
% Parameters for gas-phase calculations
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

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%---------------------------------------------------------
% Variables for calculation of q_gas according to Batstone
global P_gas
global p_gas_h2
global p_gas_ch4
global p_gas_co2
global q_gas

%/-----------------------------------------------/
%/ CALCULATIONS SECTION /
%/-----------------------------------------------/
%--------------------------------------------
%CALCULATION WITHOUT ANY ADJUSTMENT FOR K_H_I
%--------------------------------------------
factor = ((1.0/T_base) - (1.0/T_op))*(1/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); %T adjustment for K_a_co2
K_a_IN = (10^(-pK_a_IN_base))*exp(51965.0*factor); % T adjustment for K_a_IN
K_w = (10^(-pK_w_base))*exp(55900.0*factor); % T adjustment for K_w
K_H_h2 = K_H_h2_base*exp(-14240.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(51965.0*factor);

%reações ácidas
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)) ); % This equation is orriginate from (*) reference
r_A_11 =( k_A_BIN*(y(32)*(K_a_IN + y(36)) - K_a_IN*y(11)) ); % Note: S_nh4_ion = S_IN - S_nh3, S_nh4_ion is not S_IN

%equações de transferência de gás
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 )) );

%--------------- DONE -----------------
%------------------------
% Stoich (i) calculations
%------------------------
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;

%--------------- DONE -----------------
%------------------------
% Inhibition calculations
%------------------------
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 e Itepro dependem da concentração da presença de elementos traçoes
% se TEM o valor das constantes será 1 ou 0
% Sem TE
% se TAN <5g/L Itec4=Itepro=1
% se TAN >5g/L Itec4=Itepro=0
% Com TE
% Com TE
% se TAN <8g/L Itec4=Itepro=1
% se TAn >8g/L Itec4=Itepro=0

Itec4 = 1; %inibição por elementros traços de va e bu uptake
Itepro = 1; %inibição por elementos traços de pro uptake
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;

%--------------- DONE -----------------
%-----------------------------------------------
% Calculate reaction rates ro(1-19) of processes
%-----------------------------------------------
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);

%--------------- DONE -----------------
%----------------------
% gas flow calculations
%----------------------
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 adjustement for water vapour saturation pressure
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

%--------------- DONE -----------------
%----------------------
% Differential equations
%----------------------
% S_Su
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))));
% S_IN
y1(11) = (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_bac-N_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro18b+ro19)+(N_xc-f_xI_xc*N_I-f_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);
% S_h_ion
%---------------------------------------------------------
% Calculation of dS_H+ in the Thamsiriroj and Murphy, 2011
%---------------------------------------------------------
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_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_h2*N_bac*ro12+(N_bac-N_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro19)+(N_xc-f_xI_xc*N_I-f_sI_xc*N_I-f_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_bac-N_xc)*(ro13+ro14+ro15+ro16+ro17+ro18+ro18b+ro19)+(N_xc-f_xI_xc*N_I-f_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; % dS_H+ / dt

% X_ac2 decay of ac oxidisers
y1(37) = (q/V_liq)*(X_ac2f - y(37)) + Y_ac2*ro11b - ro18b;

clear q_gas

ad

% ---------------------------------------------------------------------|
% This file contains codes for solving DEs and examples of output graphs |
% ---------------------------------------------------------------------|
function ad(tspan,u)
format long
[t,y]=ode15s('ADDE',tspan,u);

%^^^^^^^^^^^^^^a^^^^^^^^^^^^^^^
%-----------------------------
% Global Varialbes
global Y_ac2;
global k_m_ac;
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 inhib89;
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;

%maxx - 10
%tspan=[0,maxx]

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); %ao mudar 9 para 10
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); %mudar 10 para 11
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('% in Biogas)','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');