Setting min value of variables within a function

Hi,

Please, how can i set a minimum within a function to zero?

I realise matlab is calculation with negative values and in the process returning some negative figures, so i would like to fix the min values of the variables to zero. this is a function of many variables that i am trying to solve with ode45.

function fval=MBBRFun4(t,y)
%Define the three variables
Xaob=y(1);
Xnb=y(2);
Xnsp=y(3);
Xcmx=y(4);
Xamx=y(5);
Xhan=y(6);
Xhaer=y(7);
Xs=y(8);
Sse=y(9);
Sno3=y(10);
Sno2=y(11);
Snh4=y(12);
So2=y(13);

% main functions
fval(1,1)=(0-V*Xaob/100)/V + u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - baob*(So2/(Ko2aob+So2))*Xaob - baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob;

fval(2,1)=(0-V*Xnb/100)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnb;

fval(3,1)=(0-V*Xnsp/100)/V + u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - bnsp*(So2/(Ko2nsp+So2))*Xnsp - bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp;

fval(4,1)=(0-V*Xcmx/100)/V + u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) + u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - bcmx*(So2/(Ko2cmx+So2))*Xcmx - bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx;

fval(5,1)=(0-V*Xamx/100)/V + u5*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx-bamx*(So2/(Ko2amx+So2))*Xamx - bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx;

fval(6,1)=(0-V*Xhan/100)/V + (u6*nhan*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u7*nhan*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan;

fval(7,1)=(0-V*Xhaer/100)/V + u8*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer - bhaer*(So2/(Ko2haer+So2))*Xhaer - bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3haer+Sno2+Sno3)*Xhaer;

fval(8,1)=(0-V*Xs*100)/V - KH*(Xs/(Xhan+Xhaer))/(KX+(Xs/(Xhan+Xhaer)))*(Xhan+Xhaer) + (1-fI)*(baob*Xaob+bnb*Xnb+bnsp*Xnsp+bcmx*Xcmx+bamx*Xamx+bhan*Xhan+bhaer*Xhaer);

fval(9,1)=(Qo*Sseo-Qe*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - (u7*nhan*1/YNo3han*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - (u8/Yhaer*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer) + KH*(Xs/(Xhan+Xhaer))/(KX+(Xs/(Xhan+Xhaer)))*(Xhan+Xhaer);

fval(10,1)=(Qo*Sno3o-Qe*Sno3)/V - (u6*nhan*(1-YNo3han)/(2.86*YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan + u5/1.14*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx + (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) + (u3/Ynsp*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) + (u4/Ycmx*(Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2))*Xcmx) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3)*Xnb - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp - ((1-fI)/2.86)*bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx - ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx - ((1-fI)/2.86)*bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhaer - ((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan;

fval(11,1)=(Qo*Sno2o-Qe*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan) - (((u5/Yamx)+u5/1.14)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx) + (u1/Yaob*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob) - (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) - (u3/Ynsp*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnsp) - (u4/Ycmx*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xcmx);
 
fval(12,1)=(Qo*Snh4o-Qe*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob) - (u4*(INBM+1/Ycmx)*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4))*Xcmx) - (u5*(INBM+1/Yamx)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx) - (u2*INBM*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) -(u3*INBM*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) - INBM*u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan - u7*INBM*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))*Xhan - u8*INBM*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer + (INBM-(fI*INXI))*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob + (INBM-(fI*INXI))*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnb + (INBM-(fI*INXI))*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp + (INBM-(fI*INXI))*bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx + (INBM-(fI*INXI))*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx + (INBM-(fI*INXI))*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan + (INBM-(fI*INXI))*bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3haer+Sno2+Sno3)*Xhaer+(INBM-(fI*INXI))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb + (INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bcmx*(So2/(Ko2cmx+So2))*Xcmx + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan + (INBM-(fI*INXI))*bhaer*(So2/(Ko2haer+So2))*Xhaer;

fval(13,1)=(Qo*So2o-Qe*So2)/V + (KaL*(So2G-So2))-((1-Yhaer)/Yhaer)*u8*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer-((3.43-Yaob)/Yaob)*u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - ((1.14-Ynb)/Ynb)*u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - ((1.14-Ynsp)/Ynsp)*u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - ((1.14-Ycmx)/Ycmx)*u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) - ((4.57-Ycmx)/Ycmx)*u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - (1-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb - (1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bcmx*(So2/(Ko2cmx+So2))*Xcmx - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan - (1-fI)*bhaer*(So2/(Ko2haer+So2))*Xhaer;

script

%initial conditions
y0=[0.02352;0.00048;0.0288;0.00432;0.216;0.1104;0.1104;2;2;2;55;50;0.5];
h=0.0006944444;
tSpan=[0 535];
[tSol, ySol]=ode15s(@(t, y) MBBRFun4(t,y), tSpan, y0);

Please people assist

Thanks

Answer 1

[Fabio Freschi]

Maybe I don't understand correctly: if you want to set all negative entries of fval to 0 you can add

fval(fval < 0) = 0;

at the end of your function. Otherwise follow Walter answer.

[Kosgey Kip]

Thanks a lot @Fabio and @Walter.

fval(fval < 0) = 0; should suffice for now.

Thanks again!!

[Kosgey Kip]

@Fabio, fval(fval < 0) = 0; unfortunately is not working, the values i am getting are unrealistically high!

I am therefore trying to create an events_function though i am getting errors

[Stephen23]

Fewer operations:

fval = max(0,fval);

[Kosgey Kip]

Thanks but this is also giving some ambiguously high values

Answer 2

[Walter Roberson]

https://www.mathworks.com/help/matlab/ref/odeset.html#d117e993600

You can set the NonNegative option to the list of variable indices that must not be negative.

However, the intention for this is strictly for variables that might become negative due to arithmetic round-off, and which would not give negative in infinite precision calculations. If that is not your situation, then instead of using NonNegative you should be using event functions to detect the component going negative and restarting the calculation as if you had "bounced" off of the negative section; see the ballode example.

[Kosgey Kip]

Thank you !!

[Kosgey Kip]

@Walter, please advice where i have gone wrong here, i have created an event_function besides the main function, and my script now looks like this:

y0=[0.02352;0.00048;0.0288;0.00432;0.216;0.1104;0.1104;2;2;2;55;50;0.5];

tSpan=[0 535];
options = odeset('Events',@eventFun);
[tSol, ySol]=ode15s(@(t, y) MBBRFun4(t,y), tSpan, y0, options);

events_function

function [value,isterminal,direction] = eventFun(t,y)
%Define the three variables
Xaob=y(1);
Xnb=y(2);
Xnsp=y(3);
Xcmx=y(4);
Xamx=y(5);
Xhan=y(6);
Xhaer=y(7);
Xs=y(8);
Sse=y(9);
Sno3=y(10);
Sno2=y(11);
Snh4=y(12);
So2=y(13);
fval(1,1)=(0-V*Xaob/100)/V + u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - baob*(So2/(Ko2aob+So2))*Xaob - baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob;
fval(2,1)=(0-V*Xnb/100)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnb;

fval(3,1)=(0-V*Xnsp/100)/V + u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - bnsp*(So2/(Ko2nsp+So2))*Xnsp - bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp;
fval(4,1)=(0-V*Xnsp/100)/V+ u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) + u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - bcmx*(So2/(Ko2cmx+So2))*Xcmx - bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx;
fval(5,1)= (0-V*Xcmx/100)/V + u5*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx-bamx*(So2/(Ko2amx+So2))*Xamx - bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx;
fval(6,1)=(0-V*Xhan/100)/V + (u6*nhan*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u7*nhan*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan;
fval(7,1)=(0-V*Xhaer/100)/V + u8*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer - bhaer*(So2/(Ko2haer+So2))*Xhaer - bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3haer+Sno2+Sno3)*Xhaer;
fval(8,1)=(0-V*Xs)/V - KH*(Xs/(Xhan+Xhaer))/(KX+(Xs/(Xhan+Xhaer)))*(Xhan+Xhaer) + (1-fI)*(baob*Xaob+bnb*Xnb+bnsp*Xnsp+bcmx*Xcmx+bamx*Xamx+bhan*Xhan+bhaer*Xhaer);
fval(9,1)=(Qo*Sseo-Qe*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - (u7*nhan*1/YNo3han*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - (u8/Yhaer*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer) + KH*(Xs/(Xhan+Xhaer))/(KX+(Xs/(Xhan+Xhaer)))*(Xhan+Xhaer);
fval(10,1)=(Qo*Sno3o-Qe*Sno3)/V - (u6*nhan*(1-YNo3han)/(2.86*YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan + u5/1.14*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx + (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) + (u3/Ynsp*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) + (u4/Ycmx*(Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2))*Xcmx) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3)*Xnb - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp - ((1-fI)/2.86)*bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx - ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx - ((1-fI)/2.86)*bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhaer - ((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan;
fval(11,1)=(Qo*Sno2o-Qe*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan) - (((u5/Yamx)+u5/1.14)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx) + (u1/Yaob*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob) - (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) - (u3/Ynsp*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnsp) - (u4/Ycmx*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xcmx);
 fval(12,1)=(Qo*Snh4o-Qe*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob) - (u4*(INBM+1/Ycmx)*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4))*Xcmx) - (u5*(INBM+1/Yamx)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx) - (u2*INBM*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) -(u3*INBM*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) - INBM*u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan - u7*INBM*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))*Xhan - u8*INBM*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer + (INBM-(fI*INXI))*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob + (INBM-(fI*INXI))*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnb + (INBM-(fI*INXI))*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp + (INBM-(fI*INXI))*bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx + (INBM-(fI*INXI))*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx + (INBM-(fI*INXI))*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan + (INBM-(fI*INXI))*bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3haer+Sno2+Sno3)*Xhaer+(INBM-(fI*INXI))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb + (INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bcmx*(So2/(Ko2cmx+So2))*Xcmx + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan + (INBM-(fI*INXI))*bhaer*(So2/(Ko2haer+So2))*Xhaer;
fval(13,1)=(Qo*So2o-Qe*So2)/V + (KaL*(So2G-So2))-((1-Yhaer)/Yhaer)*u8*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer-((3.43-Yaob)/Yaob)*u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - ((1.14-Ynb)/Ynb)*u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - ((1.14-Ynsp)/Ynsp)*u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - ((1.14-Ycmx)/Ycmx)*u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) - ((4.57-Ycmx)/Ycmx)*u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - (1-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb - (1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bcmx*(So2/(Ko2cmx+So2))*Xcmx - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan - (1-fI)*bhaer*(So2/(Ko2haer+So2))*Xhaer;
value = y(:); 
isterminal = 13; 
direction = 0;
end

However, this throws out the following error:

Attempted to access isterminal(5); index out of bounds because numel(isterminal)=1.

Error in odezero (line 142) if any(isterminal(indzc))

Error in ode15s (line 815)

[te,ye,ie,valt,stop] =

odezero(@ntrp15s,eventFcn,eventArgs,valt,...

Error in MBBR2 (line 7)

[tSol, ySol]=ode15s(@(t, y) MBBRFun4(t,y),

tSpan, y0, options);

Please advce

[Walter Roberson]

Do not set isterminal to the index of an event that is terminal. Instead set it to a vector that indicates whether each value is terminal.

Note that for your purposes, you probably only care about terminal events, so you probably do not even need to compute the other ones.