I am struggling to set a maximum value for a output for a function. I have nested my function in another but it is giving an error.

I have the below script and main function. After nesting so as to introduce a limit for maxSnh4, it is giving an error. Any hint on how to restructure it would be highly appreciated.

%Script:

tSpan=[0 300];
y0=[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) SBR(t,y), tSpan, y0);
%plotting
figure (1)
plot(tSol,ySol(:,7),'-r'); hold on %Sse
plot(tSol,ySol(:,8),'-g'); hold on %Sno3
plot(tSol,ySol(:,9),'--r'); hold on %Sno2
plot(tSol,ySol(:,10),'-k'); hold on %nh4
plot(tSol,ySol(:,11),'-b'); hold on %So2

figure (2)
plot(tSol,ySol(:,1),'-r'); hold on %Xaob
plot(tSol,ySol(:,2),'-g'); hold on %Xnb
plot(tSol,ySol(:,3),'--r'); hold on %Xnsp
plot(tSol,ySol(:,4),'-k'); hold on %amx
plot(tSol,ySol(:,5),'-b'); hold on %Xhan
plot(tSol,ySol(:,6),'--g'); hold off %Xs

xlabel('Time (days)'); ylabel('Concentration (g/m^3)');

%Main function

% SBR simulation in MATLAB
function dy=SBR(t,y)
maxSnh4=200;
    function fval=SBR1(t,y)
        % Initialize arrays
        Xaob=y(1);
        Xnb = y(2); 
        Xnsp= y(3); 
        Xamx= y(4);
        Xhan= y(5);
        Xs= y(6);
        Sse= y(7);
        Sno3= y(8);
        Sno2= y(9);
        Snh4= y(10);
        So2= y(11);
        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;
        Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
        Vo=V-0.3*V; % volume in the reactor at the start of refilling
        % SBR operation cycle
        refilling_duration = 1.59;  % Duration of refilling phase (min)
        aeration_duration = 15;   % Duration of aeration phase (min)
        mixing_duration = 45;     % Duration of mixing phase (min)
        settling_duration = 30;   % Duration of settling phase (min)
        withdrawal_duration = 0.58;  % Duration of withdrawal phase (min)
        cycle_duration=refilling_duration+aeration_duration+mixing_duration+settling_duration+withdrawal_duration;
        t0=0+cycle_duration;
        t1=t0+15;
        t2=t1+15;
        t3=t2+45;
        t4=t3+60;
        t5=t4+15;
        %refilling
        while t>t0 && t<=t1
            % input parameters
            Sseo=300; Sno3o=0.2; Qo=220; Xso=250;
            Snh4o=650; Sno2o=0.2; So2o=0; KaL=0;
            Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
            break
        end
        
        %aeration
        while t>t1 && t<=t2
            % input parameters
            Sseo=0; Sno3o=0; Qo=0; Xso=0; Snh4o=0; Sno2o=0; So2o=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
            if Snh4<10 || Sno2>10 || So2>0.8
                KaL=0;
            elseif Snh4>10 && So2<=0.5 
                KaL=1920;
            elseif Sno2<10 && So2<0.5
                KaL=1920;
            end
            break
        end
        %mixing & settling
        while t>t2 && t<=t4
            % input parameters
            Sseo=0; Sno3o=0; Qo=0; Xso=0; Snh4o=0; Sno2o=0.2; So2o=0; KaL=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
            break
        end
        %withdrawal
        while t>t4 && t<=t5
            % input parameters
            Sseo=0; Sno3o=0; Qo=600; Xso=0; Snh4o=0; Sno2o=0.2; So2o=0; KaL=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
            break
        end
        
        fval(1)=((Vo*Xaobo-Vo*Xaob/SRT)/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)=((Vo*Xnbo-Vo*Xnb/SRT)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*((Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3))*Xnb);
        fval(3)=((Vo*Xnspo-Vo*Xnsp/SRT)/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)=((Vo*Xamxo-Vo*Xamx/SRT)/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(5)=((Vo*Xhano-Vo*Xhan/SRT)/V + (u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u6*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))+u6*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan)- bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*((Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3))*Xhan);
        fval(6)=((Qo*Xso-Qo*Xs)/V - KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan);
        
        fval(7)=((Qo*Sseo-Qo*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan -(u6*nhan*(1/YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - ((u6/Yhaer)*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2))*Xhan) + KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan);
        fval(8)=((Qo*Sno3o-Qo*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) - ((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)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan);
        fval(9)=((Qo*Sno2o-Qo*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(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/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp);
        
        fval(10)=((Qo*Snh4o-Qo*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob -(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*nhan*((Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*nhan*((Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan + (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))*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))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb +(INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan);
        fval(11)=((Qo*So2o-Qo*So2)/V + KaL*(8-So2) -(((1-Yhaer)/Yhaer)*u6*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan -(((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-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb -(1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan);
        
        if fval(10)>maxSnh4
            dy=maxSnh4;
        else
            dy=fval(10);
        end
        
    end
end

[Jan]

Please post a copy of the complete error message. Usually these messages are clear and offer a solution already.

As far as I can see, maxSnh4 is not defined before used in:

if fval(10)>maxSnh4

[KIPROTICH KOSGEY]

maxSnh4 is defined in row 3.

It is my first time coding nested functions and I am struggling to structure it properly hence the errors.

the error is: Output argument "dy" (and possibly others) not assigned a value in the execution with "SBR" function.

I just need to set fval(10) to <=200. Any hint would be greatly appreciated.

[Stephen23]

"Output argument "dy" (and possibly others) not assigned a value in the execution with "SBR" function."

Because you never define the variable DY inside of the function SBR.

You also never call the nested function SBR1, so you might as well just delete it. Your code is exactly equivalent to this:

function dy=SBR(t,y)
maxSnh4=200;
end

which clearly does not define DY anywhere.

Answer 1

[Torsten]

Instead of programming if-statements within the function routine and setting parameters depending on the actual time t, you should integrate up to t0, return control to the calling program and save results, reset parameters, integrate from t0 to t1, return control to the calling program and save results, reset parameters, integrate ... This is the usual procedure to avoid discontinuities in the ODE equations which most probably makes ode15s stop with an error message at time t ~ 107.

%Script:

tSpan=[0 300];
y0=[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) SBR(t,y), tSpan, y0);

Warning: Failure at t=1.071763e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.273737e-13) at time t.

%plotting
figure (1)
plot(tSol,ySol(:,7),'-r'); hold on %Sse
plot(tSol,ySol(:,8),'-g'); hold on %Sno3
plot(tSol,ySol(:,9),'--r'); hold on %Sno2
plot(tSol,ySol(:,10),'-k'); hold on %nh4
plot(tSol,ySol(:,11),'-b'); hold on %So2

figure (2)
plot(tSol,ySol(:,1),'-r'); hold on %Xaob
plot(tSol,ySol(:,2),'-g'); hold on %Xnb
plot(tSol,ySol(:,3),'--r'); hold on %Xnsp
plot(tSol,ySol(:,4),'-k'); hold on %amx
plot(tSol,ySol(:,5),'-b'); hold on %Xhan
plot(tSol,ySol(:,6),'--g'); hold off %Xs

xlabel('Time (days)'); ylabel('Concentration (g/m^3)');

%Main function

% SBR simulation in MATLAB
function dy=SBR(t,y)
maxSnh4=200;
    %function fval=SBR1(t,y)
        % Initialize arrays
        Xaob=y(1);
        Xnb = y(2); 
        Xnsp= y(3); 
        Xamx= y(4);
        Xhan= y(5);
        Xs= y(6);
        Sse= y(7);
        Sno3= y(8);
        Sno2= y(9);
        Snh4= y(10);
        So2= y(11);
        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;
        Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
        Vo=V-0.3*V; % volume in the reactor at the start of refilling
        % SBR operation cycle
        refilling_duration = 1.59;  % Duration of refilling phase (min)
        aeration_duration = 15;   % Duration of aeration phase (min)
        mixing_duration = 45;     % Duration of mixing phase (min)
        settling_duration = 30;   % Duration of settling phase (min)
        withdrawal_duration = 0.58;  % Duration of withdrawal phase (min)
        cycle_duration=refilling_duration+aeration_duration+mixing_duration+settling_duration+withdrawal_duration;
        t0=0+cycle_duration;
        t1=t0+15;
        t2=t1+15;
        t3=t2+45;
        t4=t3+60;
        t5=t4+15;
        %refilling
        %while t>t0 && t<=t1
        if t>t0 && t<=t1
            % input parameters
            Sseo=300; Sno3o=0.2; Qo=220; Xso=250;
            Snh4o=650; Sno2o=0.2; So2o=0; KaL=0;
            Xaobo=2; Xnbo=2; Xnspo=2; Xhano=2; Xamxo=2;
            %break
        end
        
        %aeration
        %while t>t1 && t<=t2
        if t>t1 && t<=t2
            % input parameters
            Sseo=0; Sno3o=0; Qo=0; Xso=0; Snh4o=0; Sno2o=0; So2o=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
            if Snh4<10 || Sno2>10 || So2>0.8
                KaL=0;
            elseif Snh4>10 && So2<=0.5 
                KaL=1920;
            elseif Sno2<10 && So2<0.5
                KaL=1920;
            end
            %break
        end
        %mixing & settling
        %while t>t2 && t<=t4
        if t>t2 && t<=t4
            % input parameters
            Sseo=0; Sno3o=0; Qo=0; Xso=0; Snh4o=0; Sno2o=0.2; So2o=0; KaL=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
            %break
        end
        %withdrawal
        %while t>t4 && t<=t5
        if t>t4 && t<=t5
            % input parameters
            Sseo=0; Sno3o=0; Qo=600; Xso=0; Snh4o=0; Sno2o=0.2; So2o=0; KaL=0; Xaobo=0; Xnbo=0; Xnspo=0; Xhano=0; Xamxo=0;
            %break
        end
        
        fval(1)=((Vo*Xaobo-Vo*Xaob/SRT)/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)=((Vo*Xnbo-Vo*Xnb/SRT)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*((Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3))*Xnb);
        fval(3)=((Vo*Xnspo-Vo*Xnsp/SRT)/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)=((Vo*Xamxo-Vo*Xamx/SRT)/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(5)=((Vo*Xhano-Vo*Xhan/SRT)/V + (u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u6*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))+u6*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan)- bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*((Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3))*Xhan);
        fval(6)=((Qo*Xso-Qo*Xs)/V - KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan);
        
        fval(7)=((Qo*Sseo-Qo*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan -(u6*nhan*(1/YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - ((u6/Yhaer)*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2))*Xhan) + KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan);
        fval(8)=((Qo*Sno3o-Qo*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) - ((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)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan);
        fval(9)=((Qo*Sno2o-Qo*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(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/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp);
        
        fval(10)=((Qo*Snh4o-Qo*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob -(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*nhan*((Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*nhan*((Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan + (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))*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))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb +(INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan);
        fval(11)=((Qo*So2o-Qo*So2)/V + KaL*(8-So2) -(((1-Yhaer)/Yhaer)*u6*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan -(((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-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb -(1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan);
        
        if fval(10)>maxSnh4
            fval(10)=maxSnh4;
        %else
        %    dy=fval(10);
        end
        dy = fval.';
    end
%end

[KIPROTICH KOSGEY]

Thank you for the insights. I will re-structure it further. However, I have coded for fval(10)<=200 but it is still going over this value as is evident in figure(1).

Please advise on this on this.

[Torsten]

The figure doesn't display "fval", but "y".

And the mass matrix must be of size 11x11, not 55x55.

But since you don't use the options settings, MATLAB does not throw an error.