Index exceeds the number of array elements. Index must not exceed 1

I keep gettting the above error despite making adjustments. Please can someone help?

%script
tSpan=[0 5];
y0=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
M1=eye(55);
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'Mass',M1,'MassSingular', 'maybe','MStateDependence','strong','InitialStep',1e-10, 'MStateDependence','strong');
[tSol, ySol]=ode15s(@(t,y) new(t,y), tSpan, y0, options);

Index exceeds the number of array elements. Index must not exceed 1.

Error in solution>new (line 131)
Xaob1(ii)=(((Vo*Xaobo-Vo*Xaob5(ii))/SRT)/V + u1*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii)))*Xaob1(ii) - baob*(So21(ii)/(Ko2aob+So21(ii)))*Xaob1(ii) - baob*naob*((Ko2aob/(Ko2aob+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3aob+Sno21(ii)+Sno31(ii)))*Xaob1(ii));

Error in solution>@(t,y)new(t,y) (line 6)
[tSol, ySol]=ode15s(@(t,y) new(t,y), tSpan, y0, options);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

%plotting
figure (1)
plot(tSol,ySol(:,51),'-r'); hold on %Sse
plot(tSol,ySol(:,52),'-g'); hold on %Sno3
plot(tSol,ySol(:,53),'--r'); hold on %Sno2
plot(tSol,ySol(:,54),'-k'); hold on %nh4
plot(tSol,ySol(:,55),'-b'); hold on %So2
figure (2)
plot(tSol,ySol(:,45),'-r'); hold on %Xaob
plot(tSol,ySol(:,46),'-g'); hold on %Xnb
plot(tSol,ySol(:,47),'--r'); hold on %Xnsp
plot(tSol,ySol(:,48),'-k'); hold on %amx
plot(tSol,ySol(:,49),'-b'); hold on %Xhan
plot(tSol,ySol(:,50),'--g'); hold off %Xs
xlabel('Time (days)'); ylabel('Concentration (g/m^3)');

%Function
% SBR simulation in MATLAB
function fval=new(t,y)
% Initialize arrays
Xaob1 = y(1); 
Xnb1 = y(2); 
Xnsp1= y(3); 
Xamx1= y(4);
Xhan1= y(5);
Xs1= y(6);
Sse1= y(7);
Sno31= y(8);
Sno21= y(9);
Snh41= y(10);
So21= y(11);
Xaob2 = y(12); 
Xnb2 = y(13); 
Xnsp2= y(14); 
Xamx2= y(15);
Xhan2= y(16);
Xs2= y(17);
Sse2= y(18);
Sno32= y(19);
Sno22= y(20);
Snh42= y(21);
So22= y(22);
Xaob3 = y(23); 
Xnb3 = y(24); 
Xnsp3= y(25); 
Xamx3= y(26);
Xhan3= y(27);
Xs3= y(28);
Sse3= y(29);
Sno33= y(30);
Sno23= y(31);
Snh43= y(32);
So23= y(33);
Xaob4 = y(34); 
Xnb4 = y(35); 
Xnsp4= y(36); 
Xamx4= y(37);
Xhan4= y(38);
Xs4= y(39);
Sse4= y(40);
Sno34= y(41);
Sno24= y(42);
Snh44= y(43);
So24= y(44);
Xaob5 = y(45); 
Xnb5 = y(46); 
Xnsp5= y(47); 
Xamx5= y(48);
Xhan5= y(49);
Xs5= y(50);
Sse5= y(51);
Sno35= y(52);
Sno25= y(53);
Snh45= y(54);
So25= y(55);
T=303;
%aob
O1=68/(0.00831*293*(273+T)); u1=0.054*exp(O1*(T-293)); Ko2aob=0.6; Knh4aob=2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.047*exp(O2*(T-293)); Ko2nb=2.2; Kno2nb=0.5*1.375; Ynb=0.041;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.69*exp(O3*(T-293)); Ko2nsp=0.33; Kno2nsp=0.52; Ynsp=0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=0.0022*exp(O5*(T-293)); Kno2amx=0.05; Knh4amx=0.07; Yamx=0.159; Ko2amx=0.01; Kno3amx=0.5;
%han 
O6=0.069;u6=0.5*7.2*exp((T-293)*0.069); KSsehan=2; Kno2han=0.5;Kno3han=0.5; YNo2han=0.3; Ko2han=0.2;  
YNo3han=0.43; Yhaer=0.54; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=0.0025*exp(O6*(T-293));bnsp=bnb; bamx=0.00013*exp(O6*(T-293)); bhan=0.008*exp(O6*(T-293));
INBM=0.07; fI=0.1; namx=0.5; nhan=0.6;naob=0.5; nnb=0.5; nnsp=0.5;
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Qo=116.67; Xso=250;
Snh4o=650; Sno2o=0.2; So2o=0;
%oxygen transfer parameters
KaL=1920;
V=1400; SRT=100;
% Time
t_end = 5 * 24 * 60;     % Total simulation time (min)
dt = 0.1;        % Time step (min)
t = 0:dt:t_end;  % Time vector
N=(t_end-0)/dt;
Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
Vo=V-0.3*V; % volume in the reactor at the start of refilling
t1=Vo/Qo;
% SBR operation cycle
refilling_duration = t1;  % Duration of refilling phase (min)
aeration_duration = 15;   % Duration of aeration phase (min)
mixing_duration = 45;     % Duration of mixing phase (min)
settling_duration = 60;   % Duration of settling phase (min)
withdrawal_duration = t1;  % Duration of withdrawal phase (min)
% Calculate the total duration of each cycle
cycle_duration = refilling_duration + aeration_duration + mixing_duration + settling_duration + withdrawal_duration;
for i = 1:N
    % Determine the current cycle
    current_cycle = ceil(t(i)/cycle_duration);
    
    % Determine the phase of the cycle
    cycle_phase = mod(t(i), cycle_duration);
    %index
    ii=i+1;
    % Refilling phase
    if cycle_phase >= 0 && cycle_phase < refilling_duration
        % Reactor dynamics
        Xaob1(ii)=(((Vo*Xaobo-Vo*Xaob5(ii))/SRT)/V + u1*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii)))*Xaob1(ii) - baob*(So21(ii)/(Ko2aob+So21(ii)))*Xaob1(ii) - baob*naob*((Ko2aob/(Ko2aob+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3aob+Sno21(ii)+Sno31(ii)))*Xaob1(ii));
        Xnb1(ii)=(((Vo*Xnbo-Vo*Xnb5(ii))/SRT)/V + u2*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii) - bnb*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii) - bnb*nnb*((Ko2nb/(Ko2nb+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nb+Sno21(ii)+Sno31(ii)))*Xnb1(ii));
        Xnsp1(ii)=(((Vo*Xnspo-Vo*Xnsp5(ii))/SRT)/V + u3*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii) - bnsp*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii) - bnsp*nnsp*((Ko2nsp/(Ko2nsp+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nsp+Sno21(ii)+Sno31(ii)))*Xnsp1(ii));
        Xamx1(ii)=(((Vo*Xamxo-Vo*Xamx5(ii))/SRT)/V + u5*((Snh41(ii)/(Knh4amx+Snh41(ii)))*(Sno21(ii)/(Kno2amx+Sno21(ii)))*(Ko2amx/(Ko2amx+So21(ii))))*Xamx1(ii)-bamx*(So21(ii)/(Ko2amx+So21(ii)))*Xamx1(ii) - bamx*namx*((Ko2amx/(Ko2amx+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3amx+Sno21(ii)+Sno31(ii)))*Xamx1(ii));
        Xhan1(ii)=(((Vo*Xhano-Vo*Xhan5(ii))/SRT)/V + (u6*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno21(ii)/(Kno2han+Sno21(ii)))*(Ko2han/(Ko2han+So21(ii))) + u6*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno31(ii)/(Kno3han+Sno31(ii)))*(Ko2han/(Ko2han+So21(ii)))+u6*((Sse1(ii)/(KSsehan+Sse1(ii)))*(So21(ii)/(Ko2han+So21(ii))))*Xhan1(ii))- bhan*(So21(ii)/(Ko2han+So21(ii)))*Xhan1(ii) - bhan*nhan*((Ko2han/(Ko2han+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3han+Sno21(ii)+Sno31(ii)))*Xhan1(ii));
        Xs1(ii)=(((Vo*Xso-Vo*Xs5(ii))/SRT)/V - KH*((Xs1(ii)/Xhan1(ii))/(KX+(Xs1(ii)/Xhan1(ii))))*Xhan1(ii));
        
        Sse1(ii)=((Qo*Sseo-Qo*Sse5(ii))/V - (u6*nhan*(1/YNo2han)*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno21(ii)/(Kno2han+Sno21(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) -(u6*nhan*(1/YNo3han)*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno31(ii)/(Kno3han+Sno31(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) - ((u6/Yhaer)*(Sse1(ii)/(KSsehan+Sse1(ii)))*(So21(ii)/(Ko2han+So21(ii)))*Xhan1(ii)) + KH*((Xs1(ii)/Xhan1(ii))/(KX+(Xs1(ii)/Xhan1(ii))))*Xhan1(ii));
        Sno31(ii)=((Qo*Sno3o-Qo*Sno35(ii))/V - (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno31(ii)/(Kno3han+Sno31(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) + (u5/1.14)*(Snh41(ii)/(Knh4amx+Snh41(ii)))*(Sno21(ii)/(Kno2amx+Sno21(ii)))*(Ko2amx/(Ko2amx+So21(ii)))*Xamx1(ii) + (u2/Ynb*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii)) + ((u3/Ynsp)*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii)) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3aob+Sno21(ii)+Sno31(ii))*Xaob1(ii) - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nb+Sno21(ii)+Sno31(ii))*Xnb1(ii) - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nsp+Sno21(ii)+Sno31(ii))*Xnsp1(ii) -  ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3amx+Sno21(ii)+Sno31(ii))*Xamx1(ii)-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3han+Sno21(ii)+Sno31(ii))*Xhan1(ii));
        Sno21(ii)=((Qo*Sno2o-Qo*Sno25(ii))/V - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno21(ii)/(Kno2han+Sno21(ii)))*(Ko2han/(Ko2han+So21(ii)))*Xhan1(ii)) - (((u5/Yamx)+(u5/1.14))*(Snh41(ii)/(Knh4amx+Snh41(ii)))*(Sno21(ii)/(Kno2amx+Sno21(ii)))*(Ko2amx/(Ko2amx+So21(ii))))*Xamx1(ii) + ((u1/Yaob)*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii))))*Xaob1(ii) - ((u2/Ynb)*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii)) - ((u3/Ynsp)*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii))))*Xnsp1(ii));
        
        Snh41(ii)=((Qo*Snh4o-Qo*Snh45(ii))/V - (u1*(INBM+1/Yaob)*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii))))*Xaob1(ii) -(u5*(INBM+1/Yamx)*(Snh41(ii)/(Knh4amx+Snh41(ii)))*(Sno21(ii)/(Kno2amx+Sno21(ii)))*(Ko2amx/(Ko2amx+So21(ii))))*Xamx1(ii) - (u2*INBM*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii))))*Xnb1(ii) -(u3*INBM*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii))))*Xnsp1(ii) - INBM*u6*nhan*((Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno21(ii)/(Kno2han+Sno21(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) - u6*INBM*nhan*((Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno31(ii)/(Kno3han+Sno31(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) - u6*INBM*((Sse1(ii)/(KSsehan+Sse1(ii)))*(So21(ii)/(Ko2han+So21(ii))))*Xhan1(ii) + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3aob+Sno21(ii)+Sno31(ii)))*Xaob1(ii) + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nsp+Sno21(ii)+Sno31(ii)))*Xnb1(ii) +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nsp+Sno21(ii)+Sno31(ii)))*Xnsp1(ii) + (INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3amx+Sno21(ii)+Sno31(ii)))*Xamx1(ii) + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3han+Sno21(ii)+Sno31(ii)))*Xhan1(ii) +(INBM-(fI*INXI))*baob*(So21(ii)/(Ko2aob+So21(ii)))*Xaob1(ii) + (INBM-(fI*INXI))*bnb*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii) +(INBM-(fI*INXI))*bnsp*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii) + (INBM-(fI*INXI))*bamx*(So21(ii)/(Ko2amx+So21(ii)))*Xamx1(ii) + (INBM-(fI*INXI))*bhan*(So21(ii)/(Ko2han+So21(ii)))*Xhan1(ii));
        So21(ii)=((Qo*So2o-Qo*So25(ii))/V -(((1-Yhaer)/Yhaer)*u6*(Sse1(ii)/(KSsehan+Sse1(ii)))*(So21(ii)/(Ko2han+So21(ii))))*Xhan1(ii) -(((3.43-Yaob)/Yaob)*u1*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii))))*Xaob1(ii) - (((1.14-Ynb)/Ynb)*u2*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii))))*Xnb1(ii) - (((1.14-Ynsp)/Ynsp)*u3*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii))))*Xnsp1(ii) - (1-fI)*baob*(So21(ii)/(Ko2aob+So21(ii)))*Xaob1(ii) - (1-fI)*bnb*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii) -(1-fI)*bnsp*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii) - (1-fI)*bamx*(So21(ii)/(Ko2amx+So21(ii)))*Xamx1(ii) - (1-fI)*bhan*(So21(ii)/(Ko2han+So21(ii)))*Xhan1(ii));
        
    end
    
    % Aeration phase
    if cycle_phase >= refilling_duration && cycle_phase < refilling_duration + aeration_duration
        % Do aeration-related calculations or modifications here
        Xaob2(ii)=( u1*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii)))*Xaob1(ii) - baob*(So21(ii)/(Ko2aob+So21(ii)))*Xaob1(ii) - baob*naob*((Ko2aob/(Ko2aob+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3aob+Sno21(ii)+Sno31(ii)))*Xaob1(ii));
        Xnb2(ii)=( u2*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii) - bnb*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii) - bnb*nnb*((Ko2nb/(Ko2nb+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nb+Sno21(ii)+Sno31(ii)))*Xnb1(ii));
        Xnsp2(ii)=( u3*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii) - bnsp*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii) - bnsp*nnsp*((Ko2nsp/(Ko2nsp+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nsp+Sno21(ii)+Sno31(ii)))*Xnsp1(ii));
        Xamx2(ii)=( u5*((Snh41(ii)/(Knh4amx+Snh41(ii)))*(Sno21(ii)/(Kno2amx+Sno21(ii)))*(Ko2amx/(Ko2amx+So21(ii))))*Xamx1(ii)-bamx*(So21(ii)/(Ko2amx+So21(ii)))*Xamx1(ii) - bamx*namx*((Ko2amx/(Ko2amx+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3amx+Sno21(ii)+Sno31(ii)))*Xamx1(ii));
        Xhan2(ii)=( (u6*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno21(ii)/(Kno2han+Sno21(ii)))*(Ko2han/(Ko2han+So21(ii))) + u6*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno31(ii)/(Kno3han+Sno31(ii)))*(Ko2han/(Ko2han+So21(ii)))+u6*((Sse1(ii)/(KSsehan+Sse1(ii)))*(So21(ii)/(Ko2han+So21(ii))))*Xhan1(ii))- bhan*(So21(ii)/(Ko2han+So21(ii)))*Xhan1(ii) - bhan*nhan*((Ko2han/(Ko2han+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3han+Sno21(ii)+Sno31(ii)))*Xhan1(ii));
        Xs2(ii)=( - KH*((Xs1(ii)/Xhan1(ii))/(KX+(Xs1(ii)/Xhan1(ii))))*Xhan1(ii));
        
        Sse2(ii)=( - (u6*nhan*(1/YNo2han)*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno21(ii)/(Kno2han+Sno21(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) -(u6*nhan*(1/YNo3han)*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno31(ii)/(Kno3han+Sno31(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) - ((u6/Yhaer)*(Sse1(ii)/(KSsehan+Sse1(ii)))*(So21(ii)/(Ko2han+So21(ii)))*Xhan1(ii)) + KH*((Xs1(ii)/Xhan1(ii))/(KX+(Xs1(ii)/Xhan1(ii))))*Xhan1(ii));
        Sno32(ii)=( - (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno31(ii)/(Kno3han+Sno31(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) + (u5/1.14)*(Snh41(ii)/(Knh4amx+Snh41(ii)))*(Sno21(ii)/(Kno2amx+Sno21(ii)))*(Ko2amx/(Ko2amx+So21(ii)))*Xamx1(ii) + (u2/Ynb*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii)) + ((u3/Ynsp)*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii)) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3aob+Sno21(ii)+Sno31(ii))*Xaob1(ii) - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nb+Sno21(ii)+Sno31(ii))*Xnb1(ii) - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nsp+Sno21(ii)+Sno31(ii))*Xnsp1(ii) -  ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3amx+Sno21(ii)+Sno31(ii))*Xamx1(ii)-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3han+Sno21(ii)+Sno31(ii))*Xhan1(ii));
        Sno22(ii)=( - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno21(ii)/(Kno2han+Sno21(ii)))*(Ko2han/(Ko2han+So21(ii)))*Xhan1(ii)) - (((u5/Yamx)+(u5/1.14))*(Snh41(ii)/(Knh4amx+Snh41(ii)))*(Sno21(ii)/(Kno2amx+Sno21(ii)))*(Ko2amx/(Ko2amx+So21(ii))))*Xamx1(ii) + ((u1/Yaob)*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii))))*Xaob1(ii) - ((u2/Ynb)*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii)) - ((u3/Ynsp)*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii))))*Xnsp1(ii));
        
        Snh42(ii)=( - (u1*(INBM+1/Yaob)*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii))))*Xaob1(ii) -(u5*(INBM+1/Yamx)*(Snh41(ii)/(Knh4amx+Snh41(ii)))*(Sno21(ii)/(Kno2amx+Sno21(ii)))*(Ko2amx/(Ko2amx+So21(ii))))*Xamx1(ii) - (u2*INBM*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii))))*Xnb1(ii) -(u3*INBM*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii))))*Xnsp1(ii) - INBM*u6*nhan*((Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno21(ii)/(Kno2han+Sno21(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) - u6*INBM*nhan*((Sse1(ii)/(KSsehan+Sse1(ii)))*(Sno31(ii)/(Kno3han+Sno31(ii)))*(Ko2han/(Ko2han+So21(ii))))*Xhan1(ii) - u6*INBM*((Sse1(ii)/(KSsehan+Sse1(ii)))*(So21(ii)/(Ko2han+So21(ii))))*Xhan1(ii) + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3aob+Sno21(ii)+Sno31(ii)))*Xaob1(ii) + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nsp+Sno21(ii)+Sno31(ii)))*Xnb1(ii) +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3nsp+Sno21(ii)+Sno31(ii)))*Xnsp1(ii) + (INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3amx+Sno21(ii)+Sno31(ii)))*Xamx1(ii) + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So21(ii)))*(Sno21(ii)+Sno31(ii))/(Kno3han+Sno21(ii)+Sno31(ii)))*Xhan1(ii) +(INBM-(fI*INXI))*baob*(So21(ii)/(Ko2aob+So21(ii)))*Xaob1(ii) + (INBM-(fI*INXI))*bnb*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii) +(INBM-(fI*INXI))*bnsp*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii) + (INBM-(fI*INXI))*bamx*(So21(ii)/(Ko2amx+So21(ii)))*Xamx1(ii) + (INBM-(fI*INXI))*bhan*(So21(ii)/(Ko2han+So21(ii)))*Xhan1(ii));
        So22(ii)=(KaL*(8-So22(ii)) -(((1-Yhaer)/Yhaer)*u6*(Sse1(ii)/(KSsehan+Sse1(ii)))*(So21(ii)/(Ko2han+So21(ii))))*Xhan1(ii) -(((3.43-Yaob)/Yaob)*u1*(So21(ii)/(Ko2aob+So21(ii)))*(Snh41(ii)/(Knh4aob+Snh41(ii))))*Xaob1(ii) - (((1.14-Ynb)/Ynb)*u2*(Sno21(ii)/(Kno2nb+Sno21(ii)))*(So21(ii)/(Ko2nb+So21(ii))))*Xnb1(ii) - (((1.14-Ynsp)/Ynsp)*u3*(Sno21(ii)/(Kno2nsp+Sno21(ii)))*(So21(ii)/(Ko2nsp+So21(ii))))*Xnsp1(ii) - (1-fI)*baob*(So21(ii)/(Ko2aob+So21(ii)))*Xaob1(ii) - (1-fI)*bnb*(So21(ii)/(Ko2nb+So21(ii)))*Xnb1(ii) -(1-fI)*bnsp*(So21(ii)/(Ko2nsp+So21(ii)))*Xnsp1(ii) - (1-fI)*bamx*(So21(ii)/(Ko2amx+So21(ii)))*Xamx1(ii) - (1-fI)*bhan*(So21(ii)/(Ko2han+So21(ii)))*Xhan1(ii));
        
    end
    
    % Mixing phase
    if cycle_phase >= refilling_duration + aeration_duration && cycle_phase < refilling_duration + aeration_duration + mixing_duration
        % Do mixing-related calculations or modifications here
        
        Xaob3(ii)=( u1*(So22(ii)/(Ko2aob+So22(ii)))*(Snh42(ii)/(Knh4aob+Snh42(ii)))*Xaob2(ii) - baob*(So22(ii)/(Ko2aob+So22(ii)))*Xaob2(ii) - baob*naob*((Ko2aob/(Ko2aob+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3aob+Sno22(ii)+Sno32(ii)))*Xaob2(ii));
        Xnb3(ii)=( u2*(Sno22(ii)/(Kno2nb+Sno22(ii)))*(So22(ii)/(Ko2nb+So22(ii)))*Xnb2(ii) - bnb*(So22(ii)/(Ko2nb+So22(ii)))*Xnb2(ii) - bnb*nnb*((Ko2nb/(Ko2nb+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3nb+Sno22(ii)+Sno32(ii)))*Xnb2(ii));
        Xnsp3(ii)=( u3*(Sno22(ii)/(Kno2nsp+Sno22(ii)))*(So22(ii)/(Ko2nsp+So22(ii)))*Xnsp2(ii) - bnsp*(So22(ii)/(Ko2nsp+So22(ii)))*Xnsp2(ii) - bnsp*nnsp*((Ko2nsp/(Ko2nsp+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3nsp+Sno22(ii)+Sno32(ii)))*Xnsp2(ii));
        Xamx3(ii)=( u5*((Snh42(ii)/(Knh4amx+Snh42(ii)))*(Sno22(ii)/(Kno2amx+Sno22(ii)))*(Ko2amx/(Ko2amx+So22(ii))))*Xamx2(ii)-bamx*(So22(ii)/(Ko2amx+So22(ii)))*Xamx2(ii) - bamx*namx*((Ko2amx/(Ko2amx+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3amx+Sno22(ii)+Sno32(ii)))*Xamx2(ii));
        Xhan3(ii)=( (u6*(Sse2(ii)/(KSsehan+Sse2(ii)))*(Sno22(ii)/(Kno2han+Sno22(ii)))*(Ko2han/(Ko2han+So22(ii))) + u6*(Sse2(ii)/(KSsehan+Sse2(ii)))*(Sno32(ii)/(Kno3han+Sno32(ii)))*(Ko2han/(Ko2han+So22(ii)))+u6*((Sse2(ii)/(KSsehan+Sse2(ii)))*(So22(ii)/(Ko2han+So22(ii))))*Xhan2(ii))- bhan*(So22(ii)/(Ko2han+So22(ii)))*Xhan2(ii) - bhan*nhan*((Ko2han/(Ko2han+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3han+Sno22(ii)+Sno32(ii)))*Xhan2(ii));
        Xs3(ii)=( - KH*((Xs2(ii)/Xhan2(ii))/(KX+(Xs2(ii)/Xhan2(ii))))*Xhan2(ii));
        
        Sse3(ii)=( - (u6*nhan*(1/YNo2han)*(Sse2(ii)/(KSsehan+Sse2(ii)))*(Sno22(ii)/(Kno2han+Sno22(ii)))*(Ko2han/(Ko2han+So22(ii))))*Xhan2(ii) -(u6*nhan*(1/YNo3han)*(Sse2(ii)/(KSsehan+Sse2(ii)))*(Sno32(ii)/(Kno3han+Sno32(ii)))*(Ko2han/(Ko2han+So22(ii))))*Xhan2(ii) - ((u6/Yhaer)*(Sse2(ii)/(KSsehan+Sse2(ii)))*(So22(ii)/(Ko2han+So22(ii)))*Xhan2(ii)) + KH*((Xs2(ii)/Xhan2(ii))/(KX+(Xs2(ii)/Xhan2(ii))))*Xhan2(ii));
        Sno33(ii)=( - (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse2(ii)/(KSsehan+Sse2(ii)))*(Sno32(ii)/(Kno3han+Sno32(ii)))*(Ko2han/(Ko2han+So22(ii))))*Xhan2(ii) + (u5/1.14)*(Snh42(ii)/(Knh4amx+Snh42(ii)))*(Sno22(ii)/(Kno2amx+Sno22(ii)))*(Ko2amx/(Ko2amx+So22(ii)))*Xamx2(ii) + (u2/Ynb*(Sno22(ii)/(Kno2nb+Sno22(ii)))*(So22(ii)/(Ko2nb+So22(ii)))*Xnb2(ii)) + ((u3/Ynsp)*(Sno22(ii)/(Kno2nsp+Sno22(ii)))*(So22(ii)/(Ko2nsp+So22(ii)))*Xnsp2(ii)) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3aob+Sno22(ii)+Sno32(ii))*Xaob2(ii) - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3nb+Sno22(ii)+Sno32(ii))*Xnb2(ii) - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3nsp+Sno22(ii)+Sno32(ii))*Xnsp2(ii) -  ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3amx+Sno22(ii)+Sno32(ii))*Xamx2(ii)-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3han+Sno22(ii)+Sno32(ii))*Xhan2(ii));
        Sno23(ii)=( - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse2(ii)/(KSsehan+Sse2(ii)))*(Sno22(ii)/(Kno2han+Sno22(ii)))*(Ko2han/(Ko2han+So22(ii)))*Xhan2(ii)) - (((u5/Yamx)+(u5/1.14))*(Snh42(ii)/(Knh4amx+Snh42(ii)))*(Sno22(ii)/(Kno2amx+Sno22(ii)))*(Ko2amx/(Ko2amx+So22(ii))))*Xamx2(ii) + ((u1/Yaob)*(So22(ii)/(Ko2aob+So22(ii)))*(Snh42(ii)/(Knh4aob+Snh42(ii))))*Xaob2(ii) - ((u2/Ynb)*(Sno22(ii)/(Kno2nb+Sno22(ii)))*(So22(ii)/(Ko2nb+So22(ii)))*Xnb2(ii)) - ((u3/Ynsp)*(Sno22(ii)/(Kno2nsp+Sno22(ii)))*(So22(ii)/(Ko2nsp+So22(ii))))*Xnsp2(ii));
        
        Snh43(ii)=( - (u1*(INBM+1/Yaob)*(So22(ii)/(Ko2aob+So22(ii)))*(Snh42(ii)/(Knh4aob+Snh42(ii))))*Xaob2(ii) -(u5*(INBM+1/Yamx)*(Snh42(ii)/(Knh4amx+Snh42(ii)))*(Sno22(ii)/(Kno2amx+Sno22(ii)))*(Ko2amx/(Ko2amx+So22(ii))))*Xamx2(ii) - (u2*INBM*(Sno22(ii)/(Kno2nb+Sno22(ii)))*(So22(ii)/(Ko2nb+So22(ii))))*Xnb2(ii) -(u3*INBM*(Sno22(ii)/(Kno2nsp+Sno22(ii)))*(So22(ii)/(Ko2nsp+So22(ii))))*Xnsp2(ii) - INBM*u6*nhan*((Sse2(ii)/(KSsehan+Sse2(ii)))*(Sno22(ii)/(Kno2han+Sno22(ii)))*(Ko2han/(Ko2han+So22(ii))))*Xhan2(ii) - u6*INBM*nhan*((Sse2(ii)/(KSsehan+Sse2(ii)))*(Sno32(ii)/(Kno3han+Sno32(ii)))*(Ko2han/(Ko2han+So22(ii))))*Xhan2(ii) - u6*INBM*((Sse2(ii)/(KSsehan+Sse2(ii)))*(So22(ii)/(Ko2han+So22(ii))))*Xhan2(ii) + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3aob+Sno22(ii)+Sno32(ii)))*Xaob2(ii) + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3nsp+Sno22(ii)+Sno32(ii)))*Xnb2(ii) +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3nsp+Sno22(ii)+Sno32(ii)))*Xnsp2(ii) + (INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3amx+Sno22(ii)+Sno32(ii)))*Xamx2(ii) + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So22(ii)))*(Sno22(ii)+Sno32(ii))/(Kno3han+Sno22(ii)+Sno32(ii)))*Xhan2(ii) +(INBM-(fI*INXI))*baob*(So22(ii)/(Ko2aob+So22(ii)))*Xaob2(ii) + (INBM-(fI*INXI))*bnb*(So22(ii)/(Ko2nb+So22(ii)))*Xnb2(ii) +(INBM-(fI*INXI))*bnsp*(So22(ii)/(Ko2nsp+So22(ii)))*Xnsp2(ii) + (INBM-(fI*INXI))*bamx*(So22(ii)/(Ko2amx+So22(ii)))*Xamx2(ii) + (INBM-(fI*INXI))*bhan*(So22(ii)/(Ko2han+So22(ii)))*Xhan2(ii));
        So23(ii)=(-(((1-Yhaer)/Yhaer)*u6*(Sse2(ii)/(KSsehan+Sse2(ii)))*(So22(ii)/(Ko2han+So22(ii))))*Xhan2(ii) -(((3.43-Yaob)/Yaob)*u1*(So22(ii)/(Ko2aob+So22(ii)))*(Snh42(ii)/(Knh4aob+Snh42(ii))))*Xaob2(ii) - (((1.14-Ynb)/Ynb)*u2*(Sno22(ii)/(Kno2nb+Sno22(ii)))*(So22(ii)/(Ko2nb+So22(ii))))*Xnb2(ii) - (((1.14-Ynsp)/Ynsp)*u3*(Sno22(ii)/(Kno2nsp+Sno22(ii)))*(So22(ii)/(Ko2nsp+So22(ii))))*Xnsp2(ii) - (1-fI)*baob*(So22(ii)/(Ko2aob+So22(ii)))*Xaob2(ii) - (1-fI)*bnb*(So22(ii)/(Ko2nb+So22(ii)))*Xnb2(ii) -(1-fI)*bnsp*(So22(ii)/(Ko2nsp+So22(ii)))*Xnsp2(ii) - (1-fI)*bamx*(So22(ii)/(Ko2amx+So22(ii)))*Xamx2(ii) - (1-fI)*bhan*(So22(ii)/(Ko2han+So22(ii)))*Xhan2(ii));
    end
    
    % Settling phase
    if cycle_phase >= refilling_duration + aeration_duration + mixing_duration && cycle_phase < refilling_duration + aeration_duration + mixing_duration + settling_duration
        % Do settling-related calculations or modifications here
        Xaob4(ii)=( u1*(So23(ii)/(Ko2aob+So23(ii)))*(Snh43(ii)/(Knh4aob+Snh43(ii)))*Xaob3(ii) - baob*(So23(ii)/(Ko2aob+So23(ii)))*Xaob3(ii) - baob*naob*((Ko2aob/(Ko2aob+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3aob+Sno23(ii)+Sno33(ii)))*Xaob3(ii));
        Xnb4(ii)=( u2*(Sno23(ii)/(Kno2nb+Sno23(ii)))*(So23(ii)/(Ko2nb+So23(ii)))*Xnb3(ii) - bnb*(So23(ii)/(Ko2nb+So23(ii)))*Xnb3(ii) - bnb*nnb*((Ko2nb/(Ko2nb+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3nb+Sno23(ii)+Sno33(ii)))*Xnb3(ii));
        Xnsp4(ii)=( u3*(Sno23(ii)/(Kno2nsp+Sno23(ii)))*(So23(ii)/(Ko2nsp+So23(ii)))*Xnsp3(ii) - bnsp*(So23(ii)/(Ko2nsp+So23(ii)))*Xnsp3(ii) - bnsp*nnsp*((Ko2nsp/(Ko2nsp+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3nsp+Sno23(ii)+Sno33(ii)))*Xnsp3(ii));
        Xamx4(ii)=( u5*((Snh43(ii)/(Knh4amx+Snh43(ii)))*(Sno23(ii)/(Kno2amx+Sno23(ii)))*(Ko2amx/(Ko2amx+So23(ii))))*Xamx3(ii)-bamx*(So23(ii)/(Ko2amx+So23(ii)))*Xamx3(ii) - bamx*namx*((Ko2amx/(Ko2amx+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3amx+Sno23(ii)+Sno33(ii)))*Xamx3(ii));
        Xhan4(ii)=( (u6*(Sse3(ii)/(KSsehan+Sse3(ii)))*(Sno23(ii)/(Kno2han+Sno23(ii)))*(Ko2han/(Ko2han+So23(ii))) + u6*(Sse3(ii)/(KSsehan+Sse3(ii)))*(Sno33(ii)/(Kno3han+Sno33(ii)))*(Ko2han/(Ko2han+So23(ii)))+u6*((Sse3(ii)/(KSsehan+Sse3(ii)))*(So23(ii)/(Ko2han+So23(ii))))*Xhan3(ii))- bhan*(So23(ii)/(Ko2han+So23(ii)))*Xhan3(ii) - bhan*nhan*((Ko2han/(Ko2han+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3han+Sno23(ii)+Sno33(ii)))*Xhan3(ii));
        Xs4(ii)=( - KH*((Xs3(ii)/Xhan3(ii))/(KX+(Xs3(ii)/Xhan3(ii))))*Xhan3(ii));
        
        Sse4(ii)=( - (u6*nhan*(1/YNo2han)*(Sse3(ii)/(KSsehan+Sse3(ii)))*(Sno23(ii)/(Kno2han+Sno23(ii)))*(Ko2han/(Ko2han+So23(ii))))*Xhan3(ii) -(u6*nhan*(1/YNo3han)*(Sse3(ii)/(KSsehan+Sse3(ii)))*(Sno33(ii)/(Kno3han+Sno33(ii)))*(Ko2han/(Ko2han+So23(ii))))*Xhan3(ii) - ((u6/Yhaer)*(Sse3(ii)/(KSsehan+Sse3(ii)))*(So23(ii)/(Ko2han+So23(ii)))*Xhan3(ii)) + KH*((Xs3(ii)/Xhan3(ii))/(KX+(Xs3(ii)/Xhan3(ii))))*Xhan3(ii));
        Sno34(ii)=( - (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse3(ii)/(KSsehan+Sse3(ii)))*(Sno33(ii)/(Kno3han+Sno33(ii)))*(Ko2han/(Ko2han+So23(ii))))*Xhan3(ii) + (u5/1.14)*(Snh43(ii)/(Knh4amx+Snh43(ii)))*(Sno23(ii)/(Kno2amx+Sno23(ii)))*(Ko2amx/(Ko2amx+So23(ii)))*Xamx3(ii) + (u2/Ynb*(Sno23(ii)/(Kno2nb+Sno23(ii)))*(So23(ii)/(Ko2nb+So23(ii)))*Xnb3(ii)) + ((u3/Ynsp)*(Sno23(ii)/(Kno2nsp+Sno23(ii)))*(So23(ii)/(Ko2nsp+So23(ii)))*Xnsp3(ii)) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3aob+Sno23(ii)+Sno33(ii))*Xaob3(ii) - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3nb+Sno23(ii)+Sno33(ii))*Xnb3(ii) - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3nsp+Sno23(ii)+Sno33(ii))*Xnsp3(ii) -  ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3amx+Sno23(ii)+Sno33(ii))*Xamx3(ii)-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3han+Sno23(ii)+Sno33(ii))*Xhan3(ii));
        Sno24(ii)=( - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse3(ii)/(KSsehan+Sse3(ii)))*(Sno23(ii)/(Kno2han+Sno23(ii)))*(Ko2han/(Ko2han+So23(ii)))*Xhan3(ii)) - (((u5/Yamx)+(u5/1.14))*(Snh43(ii)/(Knh4amx+Snh43(ii)))*(Sno23(ii)/(Kno2amx+Sno23(ii)))*(Ko2amx/(Ko2amx+So23(ii))))*Xamx3(ii) + ((u1/Yaob)*(So23(ii)/(Ko2aob+So23(ii)))*(Snh43(ii)/(Knh4aob+Snh43(ii))))*Xaob3(ii) - ((u2/Ynb)*(Sno23(ii)/(Kno2nb+Sno23(ii)))*(So23(ii)/(Ko2nb+So23(ii)))*Xnb3(ii)) - ((u3/Ynsp)*(Sno23(ii)/(Kno2nsp+Sno23(ii)))*(So23(ii)/(Ko2nsp+So23(ii))))*Xnsp3(ii));
        
        Snh44(ii)=( - (u1*(INBM+1/Yaob)*(So23(ii)/(Ko2aob+So23(ii)))*(Snh43(ii)/(Knh4aob+Snh43(ii))))*Xaob3(ii) -(u5*(INBM+1/Yamx)*(Snh43(ii)/(Knh4amx+Snh43(ii)))*(Sno23(ii)/(Kno2amx+Sno23(ii)))*(Ko2amx/(Ko2amx+So23(ii))))*Xamx3(ii) - (u2*INBM*(Sno23(ii)/(Kno2nb+Sno23(ii)))*(So23(ii)/(Ko2nb+So23(ii))))*Xnb3(ii) -(u3*INBM*(Sno23(ii)/(Kno2nsp+Sno23(ii)))*(So23(ii)/(Ko2nsp+So23(ii))))*Xnsp3(ii) - INBM*u6*nhan*((Sse3(ii)/(KSsehan+Sse3(ii)))*(Sno23(ii)/(Kno2han+Sno23(ii)))*(Ko2han/(Ko2han+So23(ii))))*Xhan3(ii) - u6*INBM*nhan*((Sse3(ii)/(KSsehan+Sse3(ii)))*(Sno33(ii)/(Kno3han+Sno33(ii)))*(Ko2han/(Ko2han+So23(ii))))*Xhan3(ii) - u6*INBM*((Sse3(ii)/(KSsehan+Sse3(ii)))*(So23(ii)/(Ko2han+So23(ii))))*Xhan3(ii) + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3aob+Sno23(ii)+Sno33(ii)))*Xaob3(ii) + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3nsp+Sno23(ii)+Sno33(ii)))*Xnb3(ii) +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3nsp+Sno23(ii)+Sno33(ii)))*Xnsp3(ii) + (INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3amx+Sno23(ii)+Sno33(ii)))*Xamx3(ii) + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So23(ii)))*(Sno23(ii)+Sno33(ii))/(Kno3han+Sno23(ii)+Sno33(ii)))*Xhan3(ii) +(INBM-(fI*INXI))*baob*(So23(ii)/(Ko2aob+So23(ii)))*Xaob3(ii) + (INBM-(fI*INXI))*bnb*(So23(ii)/(Ko2nb+So23(ii)))*Xnb3(ii) +(INBM-(fI*INXI))*bnsp*(So23(ii)/(Ko2nsp+So23(ii)))*Xnsp3(ii) + (INBM-(fI*INXI))*bamx*(So23(ii)/(Ko2amx+So23(ii)))*Xamx3(ii) + (INBM-(fI*INXI))*bhan*(So23(ii)/(Ko2han+So23(ii)))*Xhan3(ii));
        So24(ii)=( -(((1-Yhaer)/Yhaer)*u6*(Sse3(ii)/(KSsehan+Sse3(ii)))*(So23(ii)/(Ko2han+So23(ii))))*Xhan3(ii) -(((3.43-Yaob)/Yaob)*u1*(So23(ii)/(Ko2aob+So23(ii)))*(Snh43(ii)/(Knh4aob+Snh43(ii))))*Xaob3(ii) - (((1.14-Ynb)/Ynb)*u2*(Sno23(ii)/(Kno2nb+Sno23(ii)))*(So23(ii)/(Ko2nb+So23(ii))))*Xnb3(ii) - (((1.14-Ynsp)/Ynsp)*u3*(Sno23(ii)/(Kno2nsp+Sno23(ii)))*(So23(ii)/(Ko2nsp+So23(ii))))*Xnsp3(ii) - (1-fI)*baob*(So23(ii)/(Ko2aob+So23(ii)))*Xaob3(ii) - (1-fI)*bnb*(So23(ii)/(Ko2nb+So23(ii)))*Xnb3(ii) -(1-fI)*bnsp*(So23(ii)/(Ko2nsp+So23(ii)))*Xnsp3(ii) - (1-fI)*bamx*(So23(ii)/(Ko2amx+So23(ii)))*Xamx3(ii) - (1-fI)*bhan*(So23(ii)/(Ko2han+So23(ii)))*Xhan3(ii));
        
    end
    % Withdrawal phase
    if cycle_phase >= refilling_duration + aeration_duration + mixing_duration + settling_duration && cycle_phase < refilling_duration + aeration_duration + mixing_duration + settling_duration + withdrawal_duration
        % Do withdrawal-related calculations or modifications here
        Xaob5(ii)=( u1*(So24(ii)/(Ko2aob+So24(ii)))*(Snh44(ii)/(Knh4aob+Snh44(ii)))*Xaob4(ii) - baob*(So24(ii)/(Ko2aob+So24(ii)))*Xaob4(ii) - baob*naob*((Ko2aob/(Ko2aob+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3aob+Sno24(ii)+Sno34(ii)))*Xaob4(ii));
        Xnb5(ii)=( u2*(Sno24(ii)/(Kno2nb+Sno24(ii)))*(So24(ii)/(Ko2nb+So24(ii)))*Xnb4(ii) - bnb*(So24(ii)/(Ko2nb+So24(ii)))*Xnb4(ii) - bnb*nnb*((Ko2nb/(Ko2nb+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3nb+Sno24(ii)+Sno34(ii)))*Xnb4(ii));
        Xnsp5(ii)=( u3*(Sno24(ii)/(Kno2nsp+Sno24(ii)))*(So24(ii)/(Ko2nsp+So24(ii)))*Xnsp4(ii) - bnsp*(So24(ii)/(Ko2nsp+So24(ii)))*Xnsp4(ii) - bnsp*nnsp*((Ko2nsp/(Ko2nsp+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3nsp+Sno24(ii)+Sno34(ii)))*Xnsp4(ii));
        Xamx5(ii)=( u5*((Snh44(ii)/(Knh4amx+Snh44(ii)))*(Sno24(ii)/(Kno2amx+Sno24(ii)))*(Ko2amx/(Ko2amx+So24(ii))))*Xamx4(ii)-bamx*(So24(ii)/(Ko2amx+So24(ii)))*Xamx4(ii) - bamx*namx*((Ko2amx/(Ko2amx+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3amx+Sno24(ii)+Sno34(ii)))*Xamx4(ii));
        Xhan5(ii)=( (u6*(Sse4(ii)/(KSsehan+Sse4(ii)))*(Sno24(ii)/(Kno2han+Sno24(ii)))*(Ko2han/(Ko2han+So24(ii))) + u6*(Sse4(ii)/(KSsehan+Sse4(ii)))*(Sno34(ii)/(Kno3han+Sno34(ii)))*(Ko2han/(Ko2han+So24(ii)))+u6*((Sse4(ii)/(KSsehan+Sse4(ii)))*(So24(ii)/(Ko2han+So24(ii))))*Xhan4(ii))- bhan*(So24(ii)/(Ko2han+So24(ii)))*Xhan4(ii) - bhan*nhan*((Ko2han/(Ko2han+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3han+Sno24(ii)+Sno34(ii)))*Xhan4(ii));
        Xs5(ii)=( - KH*((Xs54(ii)/Xhan4(ii))/(KX+(Xs54(ii)/Xhan4(ii))))*Xhan4(ii));
        
        Sse5(ii)=( - (u6*nhan*(1/YNo2han)*(Sse4(ii)/(KSsehan+Sse4(ii)))*(Sno24(ii)/(Kno2han+Sno24(ii)))*(Ko2han/(Ko2han+So24(ii))))*Xhan4(ii) -(u6*nhan*(1/YNo3han)*(Sse4(ii)/(KSsehan+Sse4(ii)))*(Sno34(ii)/(Kno3han+Sno34(ii)))*(Ko2han/(Ko2han+So24(ii))))*Xhan4(ii) - ((u6/Yhaer)*(Sse4(ii)/(KSsehan+Sse4(ii)))*(So24(ii)/(Ko2han+So24(ii)))*Xhan4(ii)) + KH*((Xs54(ii)/Xhan4(ii))/(KX+(Xs54(ii)/Xhan4(ii))))*Xhan4(ii));
        Sno35(ii)=( - (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse4(ii)/(KSsehan+Sse4(ii)))*(Sno34(ii)/(Kno3han+Sno34(ii)))*(Ko2han/(Ko2han+So24(ii))))*Xhan4(ii) + (u5/1.14)*(Snh44(ii)/(Knh4amx+Snh44(ii)))*(Sno24(ii)/(Kno2amx+Sno24(ii)))*(Ko2amx/(Ko2amx+So24(ii)))*Xamx4(ii) + (u2/Ynb*(Sno24(ii)/(Kno2nb+Sno24(ii)))*(So24(ii)/(Ko2nb+So24(ii)))*Xnb4(ii)) + ((u3/Ynsp)*(Sno24(ii)/(Kno2nsp+Sno24(ii)))*(So24(ii)/(Ko2nsp+So24(ii)))*Xnsp4(ii)) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3aob+Sno24(ii)+Sno34(ii))*Xaob4(ii) - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3nb+Sno24(ii)+Sno34(ii))*Xnb4(ii) - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3nsp+Sno24(ii)+Sno34(ii))*Xnsp4(ii) -  ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3amx+Sno24(ii)+Sno34(ii))*Xamx4(ii)-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3han+Sno24(ii)+Sno34(ii))*Xhan4(ii));
        Sno25(ii)=( - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse4(ii)/(KSsehan+Sse4(ii)))*(Sno24(ii)/(Kno2han+Sno24(ii)))*(Ko2han/(Ko2han+So24(ii)))*Xhan4(ii)) - (((u5/Yamx)+(u5/1.14))*(Snh44(ii)/(Knh4amx+Snh44(ii)))*(Sno24(ii)/(Kno2amx+Sno24(ii)))*(Ko2amx/(Ko2amx+So24(ii))))*Xamx4(ii) + ((u1/Yaob)*(So24(ii)/(Ko2aob+So24(ii)))*(Snh44(ii)/(Knh4aob+Snh44(ii))))*Xaob4(ii) - ((u2/Ynb)*(Sno24(ii)/(Kno2nb+Sno24(ii)))*(So24(ii)/(Ko2nb+So24(ii)))*Xnb4(ii)) - ((u3/Ynsp)*(Sno24(ii)/(Kno2nsp+Sno24(ii)))*(So24(ii)/(Ko2nsp+So24(ii))))*Xnsp4(ii));
        
        Snh45(ii)=( - (u1*(INBM+1/Yaob)*(So24(ii)/(Ko2aob+So24(ii)))*(Snh44(ii)/(Knh4aob+Snh44(ii))))*Xaob4(ii) -(u5*(INBM+1/Yamx)*(Snh44(ii)/(Knh4amx+Snh44(ii)))*(Sno24(ii)/(Kno2amx+Sno24(ii)))*(Ko2amx/(Ko2amx+So24(ii))))*Xamx4(ii) - (u2*INBM*(Sno24(ii)/(Kno2nb+Sno24(ii)))*(So24(ii)/(Ko2nb+So24(ii))))*Xnb4(ii) -(u3*INBM*(Sno24(ii)/(Kno2nsp+Sno24(ii)))*(So24(ii)/(Ko2nsp+So24(ii))))*Xnsp4(ii) - INBM*u6*nhan*((Sse4(ii)/(KSsehan+Sse4(ii)))*(Sno24(ii)/(Kno2han+Sno24(ii)))*(Ko2han/(Ko2han+So24(ii))))*Xhan4(ii) - u6*INBM*nhan*((Sse4(ii)/(KSsehan+Sse4(ii)))*(Sno34(ii)/(Kno3han+Sno34(ii)))*(Ko2han/(Ko2han+So24(ii))))*Xhan4(ii) - u6*INBM*((Sse4(ii)/(KSsehan+Sse4(ii)))*(So24(ii)/(Ko2han+So24(ii))))*Xhan4(ii) + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3aob+Sno24(ii)+Sno34(ii)))*Xaob4(ii) + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3nsp+Sno24(ii)+Sno34(ii)))*Xnb4(ii) +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3nsp+Sno24(ii)+Sno34(ii)))*Xnsp4(ii) + (INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3amx+Sno24(ii)+Sno34(ii)))*Xamx4(ii) + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So24(ii)))*(Sno24(ii)+Sno34(ii))/(Kno3han+Sno24(ii)+Sno34(ii)))*Xhan4(ii) +(INBM-(fI*INXI))*baob*(So24(ii)/(Ko2aob+So24(ii)))*Xaob4(ii) + (INBM-(fI*INXI))*bnb*(So24(ii)/(Ko2nb+So24(ii)))*Xnb4(ii) +(INBM-(fI*INXI))*bnsp*(So24(ii)/(Ko2nsp+So24(ii)))*Xnsp4(ii) + (INBM-(fI*INXI))*bamx*(So24(ii)/(Ko2amx+So24(ii)))*Xamx4(ii) + (INBM-(fI*INXI))*bhan*(So24(ii)/(Ko2han+So24(ii)))*Xhan4(ii));
        So25(ii)=( -(((1-Yhaer)/Yhaer)*u6*(Sse4(ii)/(KSsehan+Sse4(ii)))*(So24(ii)/(Ko2han+So24(ii))))*Xhan4(ii) -(((3.43-Yaob)/Yaob)*u1*(So24(ii)/(Ko2aob+So24(ii)))*(Snh44(ii)/(Knh4aob+Snh44(ii))))*Xaob4(ii) - (((1.14-Ynb)/Ynb)*u2*(Sno24(ii)/(Kno2nb+Sno24(ii)))*(So24(ii)/(Ko2nb+So24(ii))))*Xnb4(ii) - (((1.14-Ynsp)/Ynsp)*u3*(Sno24(ii)/(Kno2nsp+Sno24(ii)))*(So24(ii)/(Ko2nsp+So24(ii))))*Xnsp4(ii) - (1-fI)*baob*(So24(ii)/(Ko2aob+So24(ii)))*Xaob4(ii) - (1-fI)*bnb*(So24(ii)/(Ko2nb+So24(ii)))*Xnb4(ii) -(1-fI)*bnsp*(So24(ii)/(Ko2nsp+So24(ii)))*Xnsp4(ii) - (1-fI)*bamx*(So24(ii)/(Ko2amx+So24(ii)))*Xamx4(ii) - (1-fI)*bhan*(So24(ii)/(Ko2han+So24(ii)))*Xhan4(ii));
    end
    
    %fix output values to zero
    Sno31(Sno31<0)=0;Xaob1(Xaob1<0)=0;Xnb1(Xnb1<0)=0;Xnsp1(Xnsp1<0)=0;Xamx1(Xamx1<0)=0;Xhan1(Xhan1<0)=0;Xs1(Xs1<0)=0;
    Sse1(Sse1<0)=0;Sno21(Sno21<0)=0;Snh41(Snh41<0)=0;So21(So21<0)=0;
    %fix output values to zero
    Sno32(Sno32<0)=0;Xaob2(Xaob2<0)=0;Xnb2(Xnb2<0)=0;Xnsp2(Xnsp2<0)=0;Xamx2(Xamx2<0)=0;Xhan2(Xhan2<0)=0;Xs2(Xs2<0)=0;
    Sse2(Sse2<0)=0;Sno22(Sno22<0)=0;Snh42(Snh42<0)=0;So22(So22<0)=0;
    %fix output values to zero
    Sno33(Sno33<0)=0;Xaob3(Xaob3<0)=0;Xnb3(Xnb3<0)=0;Xnsp3(Xnsp3<0)=0;Xamx3(Xamx3<0)=0;Xhan3(Xhan3<0)=0;Xs3(Xs3<0)=0;
    Sse3(Sse3<0)=0;Sno23(Sno23<0)=0;Snh43(Snh43<0)=0;So23(So23<0)=0;
    %fix output values to zero
    Sno34(Sno34<0)=0;Xaob4(Xaob4<0)=0;Xnb4(Xnb4<0)=0;Xnsp4(Xnsp4<0)=0;Xamx4(Xamx4<0)=0;Xhan4(Xhan4<0)=0;Xs4(Xs4<0)=0;
    Sse4(Sse4<0)=0;Sno24(Sno24<0)=0;Snh44(Snh44<0)=0;So24(So24<0)=0;
    %fix output values to zero
    Sno35(Sno35<0)=0;Xaob5(Xaob5<0)=0;Xnb5(Xnb5<0)=0;Xnsp5(Xnsp5<0)=0;Xamx5(Xamx5<0)=0;Xhan5(Xhan5<0)=0;Xs5(Xs5<0)=0;
    Sse5(Sse5<0)=0;Sno25(Sno25<0)=0;Snh45(Snh45<0)=0;So25(So25<0)=0;
    fval(1)=Xaob1(ii) ; 
    fval(2)=Xnb1(ii) ; 
    fval(3)=Xnsp1(ii); 
    fval(4)=Xamx1(ii);
    fval(5)=Xhan1(ii);
    fval(6)=Xs1(ii);
    fval(7)=Sse1(ii);
    fval(8)=Sno31(ii);
    fval(9)=Sno21(ii);
    fval(10)=Snh41(ii);
    fval(11)=So21(ii);
    
    fval(12)=Xaob2(ii) ;
    fval(13)=Xnb2(ii) ; 
    fval(14)=Xnsp2(ii); 
    fval(15)=Xamx2(ii);
    fval(16)=Xhan2(ii);
    fval(17)=Xs2(ii);
    fval(18)=Sse2(ii);
    fval(19)=Sno32(ii);
    fval(20)=Sno22(ii);
    fval(21)=Snh42(ii);
    fval(22)=So22(ii);
    fval(23)=Xaob3(ii) ; 
    fval(24)=Xnb3(ii) ; 
    fval(25)=Xnsp3(ii); 
    fval(26)=Xamx3(ii);
    fval(27)=Xhan3(ii);
    fval(28)=Xs3(ii);
    fval(29)=Sse3(ii);
    fval(30)=Sno33(ii);
    fval(31)=Sno23(ii);
    fval(32)=Snh43(ii);
    fval(33)=So23(ii);
    fval(34)=Xaob4(ii) ;
    fval(35)=Xnb4(ii) ; 
    fval(36)=Xnsp4(ii); 
    fval(37)=Xamx4(ii);
    fval(38)=Xhan4(ii);
    fval(39)=Xs4(ii);
    fval(40)=Sse4(ii);
    fval(41)=Sno34(ii);
    fval(42)=Sno24(ii);
    fval(43)=Snh44(ii);
    fval(44)=So24(ii);
    fval(45)=Xaob5(ii) ; 
    fval(46)=Xnb5(ii) ; 
    fval(47)=Xnsp5(ii); 
    fval(48)=Xamx5(ii);
    fval(49)=Xhan5(ii);
    fval(50)=Xs5(ii);
    fval(51)=Sse5(ii);
    fval(52)=Sno35(ii);
    fval(53)=Sno25(ii);
    fval(54)=Snh45(ii);
    fval(55)=So25(ii);
    fval = fval(:);
    dbstop iferror
end
end

Answer 1

[Cris LaPierre]

This error means you are trying to index an element that does not exist in your variable. Here, your variables are ll scalars, but you are indexing them like they are vectors. For example

% Xaxob5 is a scalar
Xaob5 = y(45); 
% but your code indexes it like it is a vector
Xaob5(ii)

[KIPROTICH KOSGEY]

Thank you but I have hit a wall here, any hint on how I can sort it out?

Thanks again!

[Cris LaPierre]

The ode function is solving a single time point at a time. You appear to be trying to solve for all time points each time.

Rewrite your ode function to solve for a single time point. The solver inputs are the ouput values from the previous time step.

[KIPROTICH KOSGEY]

Thank you so much! After editing the codes, I managed to get the program running.