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

I keep getting the above error 'Index exceeds the number of array elements. Index must not exceed 1.' when I run the below program. What could be the problem?

Thanks in advance...

tstart = 0;  
tend = 180*24;   
nt = 14400000;   
xstart = 0;  
xend = 0.012;  
nx = 2000;
t = linspace(tstart,tend,nt);  
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y01 = 0.5;  y02 = 0.5; y03 = 0.5;y04 = 0.5; y05 = 0.5; y06 = 0.5;
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = 1;

y1=[y01; y02;y03;y04; y05;y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];   
%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
Qo=10; Snh4o=5; KaL=80;
%% execution
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'InitialStep',1e-10, 'NonNegative',1);
[tSol, ySol]=ode15s(@(t,y) fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL),t,y1, options);
%% feed parameters
function dydt=fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL)
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
muaob=u1.*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF)).*faob - baob.*(So2BF./(Ko2aob+So2BF)).*faob - baob*naob.*((Ko2aob./(Ko2aob+So2BF)).*(Sno2BF+Sno3BF)./(Kno3aob+Sno2BF+Sno3BF)).*faob;
munb=u2.*(Sno2BF./(Kno2nb+Sno2BF)).*(So2BF/(Ko2nb+So2BF)).*fnb - bnb.*(So2BF./(Ko2nb+So2BF)).*fnb - bnb*nnb.*((Ko2nb./(Ko2nb+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nb+Sno2BF+Sno3BF)).*fnb;
munsp=u3.*(Sno2BF./(Kno2nsp+Sno2BF)).*(So2BF/(Ko2nsp+So2BF)).*fnsp - bnsp.*(So2BF./(Ko2nsp+So2BF)).*fnsp - bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nsp+Sno2BF+Sno3BF)).*fnsp;
muamx=u5.*((Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF))).*famx - bamx.*(So2BF./(Ko2amx+So2BF)).*famx - bamx*namx.*((Ko2amx/(Ko2amx+So2BF)).*(Sno2BF+Sno3BF)./(Kno3amx+Sno2BF+Sno3BF)).*famx;
muhan=(u6.*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)) + u6.*(SseBF/(KSsehan+SseBF)).*(Sno3BF/(Kno3han+Sno3BF)).*(Ko2han/(Ko2han+So2BF))+u6.*((SseBF/(KSsehan+SseBF)).*(So2BF/(Ko2han+So2BF))).*fhan)- bhan.*(So2BF/(Ko2han+So2BF)).*fhan - bhan*nhan.*((Ko2han./(Ko2han+So2BF)).*(Sno2BF+Sno3BF)./(Kno3han+Sno2BF+Sno3BF)).*fhan;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L.*muaver+sigma;
VL=VR-AL.*(L + LL); 
%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); So2BFmin=So2BF(1); 
faobmin=exp((muaob(1)-muaver(1))*0.0003);fnbmin=exp((munb(1)-muaver(1))*0.0003);fnspmin=exp((munsp(1)-muaver(1))*0.0003);famxmin=exp((muamx(1)-muaver(1))*0.0003);fhanmin=exp((muhan(1)-muaver(1))*0.0003);
%BC-2
SseBFmax=SseBL; Sno3BFmax=Sno3BL; Sno2BFmax=Sno2BL; Snh4BFmax=Snh4BL;
So2BFmax=So2BL; 
SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax];
faob=[faobmin;faob];fnb=[fnbmin;fnb];fnsp=[fnspmin;fnsp];famx=[famxmin;famx];fhan=[fhanmin;fhan];
%dSseBLdt=zeros(nx,1);dSno3BLdt=zeros(nx,1);dSno2BLdt=zeros(nx,1);dSnh4BLdt=zeros(nx,1);dSo2BLdt=zeros(nx,1);dXsdt=zeros(nx,1);
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%dLdt=zeros(nx,1);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx
    hz=nx-(nx-1);
    z=0.012; zeta=z./L;
    U=L.*muaver*zeta;
    dSseBLdt=Qo*(Sseo-SseBL)/VL(nx)-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
    dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL(nx)-AL.*((DSno3/LL)-UL(nx)).*(Sno3BL-Sno3BF(nx));
    dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL(nx)-AL.*((DSno2/LL)-UL(nx)).*(Sno2BL-Sno2BF(nx));
    dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL(nx)-AL.*((DSnh4/LL)-UL(nx)).*(Snh4BL-Snh4BF(nx));
    dSo2BLdt=Qo.*(So2o-So2BL)./VL(nx) + KaL.*(8-So2BL) - AL.*((DSo2)-UL(nx)).*(So2BL-So2BF(nx));
    dXsdt=Qo.*(Xso-Xs)./VL(nx) - KH.*((Xs./fhan(nx))./(KX+(Xs/fhan(nx)))).*fhan(nx);
    dSseBFdt(ix)=DSse.*(SseBF(ix+1)-2.*SseBF(ix)+(SseBF(ix-1)))./(hz^2.*L(ix)^2) +zeta.*UL.*(SseBF(ix+1)-SseBF(ix))./L - (u6*nhan.*(1/YNo2han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) -(u6*nhan.*(1/YNo3han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - (u6/Yhaer).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSno3BFdt(ix)=DSno3.*(Sno3BF(ix+1)-2.*Sno3BF(ix)+(Sno3BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Sno3BF(ix+1)-Sno3BF(ix))./L - (u6*nhan*((1-YNo3han)/(2.86*YNo3han)).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) + (u5/1.14).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix))).*famx(ix) + (u2/Ynb*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) + ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix)) - ((1-fI)/2.86)*baob*naob.*(Ko2aob./(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix)).*faob(ix) - ((1-fI)/2.86)*bnb*nnb.*(Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nb+Sno2BF(ix)+Sno3BF(ix)).*fnb(ix) - ((1-fI)/2.86)*bnsp*nnsp.*(Ko2nsp/(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix)).*fnsp(ix) -  ((1-fI)/2.86)*bamx*namx.*(Ko2amx/(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix)).*famx(ix)-((1-fI)/2.86)*bhan*nhan.*(Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix)).*fhan(ix);
    dSno2BFdt(ix)=DSno2.*(Sno2BF(ix+1)-2.*Sno2BF(ix)+(Sno2BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Sno2BF(ix+1)-Sno2BF(ix))./L - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix))).*fhan(ix)) - (((u5/Yamx)+(u5/1.14)).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)/(Kno2amx+Sno2BF(ix))).*(Ko2amx/(Ko2amx+So2BF(ix)))).*famx(ix) + ((u1/Yaob).*(So2BF(ix)/(Ko2aob+So2(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - ((u2/Ynb).*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) - ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix);
    dSnh4BFdt(ix)=DSnh4.*(Snh4BF(ix+1)-2.*Snh4BF(ix)+(Snh4BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Snh4BF(ix+1)-Snh4BF(ix))./L - (u1*(INBM+1/Yaob).*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) -(u5*(INBM+1/Yamx).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix)))).*famx(ix) - (u2*INBM.*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*(ix) -(u3*INBM.*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - INBM*u6*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) + (INBM-(fI*INXI))*baob*naob.*((Ko2aob/(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb*nnb.*((Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx*namx.*((Ko2amx./(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan*nhan.*((Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix))).*fhan(ix) +(INBM-(fI*INXI))*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSo2BFdt(ix)=DSo2.*(So2BF(ix+1)-2.*So2BF(ix)+(So2BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(So2BF(ix+1)-So2BF(ix))./L - -(((1-Yhaer)/Yhaer)*u6.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) -(((3.43-Yaob)/Yaob)*u1.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - (((1.14-Ynb)/Ynb)*u2.*(Sno2BF(ix)/(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*fnb(ix) - (((1.14-Ynsp)/Ynsp)*u3.*(Sno2BF(ix)/(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - (1-fI)*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) - (1-fI)*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb -(1-fI)*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) - (1-fI)*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) - (1-fI)*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dfaobdt(ix)=(muaob-muaver).*faob(ix) - (U-zeta.*UL).*(faob(ix+1)-faob(ix))./hz;
    dfnbdt(ix)=(munb-muaver).*fnb(ix) - (U-zeta.*UL).*(fnb(ix+1)-fnb(ix))./hz;
    dfnspdt(ix)=(munsp-muaver).*fnsp(ix) - (U-zeta.*UL).*(fnsp(ix+1)-fnsp(ix))./hz;
    dfamxdt(ix)=(muamx-muaver).*famx(ix) - (U-zeta.*UL).*(famx(ix+1)-famx(ix))./hz;
    dfhandt(ix)=(muhan-muaver).*fhan(ix) - (U-zeta.*UL).*(fhan(ix+1)-fahan(ix))./hz;
    dLdt(ix)=L(ix).*muaver.*zeta+sigma; 
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [1;1;1;1;1;1;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end

[Dyuman Joshi]

1 - Please do not use dynamically named variables, it is not a good programming practice.

y01 = 0.5;  y02 = 0.5; y03 = 0.5;y04 = 0.5; y05 = 0.5; y06 = 0.5;

As Steve Vai says (or rather plays), For the love of god, do not do this, simply define an array

y0 = 0.5*ones(1,6)

y0 = 1×6

0.5000 0.5000 0.5000 0.5000 0.5000 0.5000

and use indexing.

2 - VL is a scalar i.e. it has only elements, and you are trying to access its nxth element, which is not possible. That's why you get the error.

%                        vvvvvv
dSseBLdt=Qo*(Sseo-SseBL)/VL(nx)-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));

[KIPROTICH KOSGEY]

I have edited but still the error persists.

New codes are:

y0 = 0.5*ones(1,6); y017=10^-6*ones(1,1);

y01=y0(1);y02=y0(2);y03=y0(3);y04=y0(4);y05=y0(5);y06=y0(6);y17=y017(1);

dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));

dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx)); dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));

dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));

dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));

dXsdt=Qo.*(Xso-Xs)./VL - KH.*((Xs./fhan(nx))./(KX+(Xs/fhan(nx)))).*fhan(nx);

[Dyuman Joshi]

It's the same thing as above, Ssebf has only 2 elements, and you are trying to access the nxth element, which is not possible.

You should check the size of all the variables before running your code.

And, this is still not good -

y01=y0(1);y02=y0(2);y03=y0(3);y04=y0(4);y05=y0(5);y06=y0(6);y17=y017(1);

Do this -

%Define y0 as a column vector
y0 = 0.5*ones(6,1);
y17 = 10^-6*ones(1,1);
%USE INDEXING!!!
y1=[y0; b(7:16,:).'; y17];   

[KIPROTICH KOSGEY]

I am struggling to code the BC for pdes.

%BC-1

SseBF(1)=SseBF(2); Sno3BF(1)=Sno3BF(2); Sno2BF(1)=Sno2BF(2); Snh4BF(1)=Snh4BF(2); So2BF(1)=So2BF(2);

faob(1)=exp((muaob(1)-muaver(1))*0.0003);fnb(1)=exp((munb(1)-muaver(1))*0.0003);fnsp(1)=exp((munsp(1)-muaver(1))*0.0003);famx(1)=exp((muamx(1)-muaver(1))*0.0003);fhan(1)=exp((muhan(1)-muaver(1))*0.0003);

%BC-2

SseBF(nx)=SseBL; Sno3BF(nx)=Sno3BL; Sno2BF(nx)=Sno2BL; Snh4BF(nx)=Snh4BL;

So2BF(nx)=So2BL;

[Torsten]

Loop from ix = 2:nx-1 and set the boundary conditions separately.

[KIPROTICH KOSGEY]

Thank you Torsten, you have been so helpful in the course of my project. I didn't quite understand what you mean by 'set the boundary conditions separately.'

I have tried introducing a for loop for BCs but it aint helping:

for ix=1:nx

%BC-1

SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); So2BFmin=So2BF(1);

faobmin=exp((muaob(1)-muaver(1))*0.0003);fnbmin=exp((munb(1)-muaver(1))*0.0003);fnspmin=exp((munsp(1)-muaver(1))*0.0003);famxmin=exp((muamx(1)-muaver(1))*0.0003);fhanmin=exp((muhan(1)-muaver(1))*0.0003);

%BC-2

SseBFmax=SseBL; Sno3BFmax=Sno3BL; Sno2BFmax=Sno2BL; Snh4BFmax=Snh4BL;

So2BFmax=So2BL;

end

SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax];

faob=[faobmin;faob];fnb=[fnbmin;fnb];fnsp=[fnspmin;fnsp];famx=[famxmin;famx];fhan=[fhanmin;fhan];

[Torsten]

The output of your function "fun" are time derivatives.

Thus if you want to define boundary conditions and if the function values in the boundary points are solution variables, you have to transfer their time derivatives to "ode15s".

If you want to prescribe a value for a variable and the value remains the same over time, you simply have to set d(variable)/dt = 0 for ix = 1 or ix = nx.

If you want to prescribe a value for a variable and the value changes over time (say like a function boundary(t)), you have to set d(variable)/dt = d(boundary)/dt for ix = 1 or ix = nx.

If you want to set d(variable)/dx = something as boundary condition, look at the test code we developed here:

https://uk.mathworks.com/matlabcentral/answers/1999518-how-do-i-fix-error-using-odearguments-returns-a-vector-of-length-0-but-the-length-of-initial#answer_1286132?s_tid=prof_contriblnk

[KIPROTICH KOSGEY]

Thank you so much.

I have nested another function to solve for the ODE for the left boundary conditions. However, I am getting an error regarding the indexing within the BC code:

function dydt=fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL)
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
dydt=@fval;
muaob=u1.*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF)).*faob - baob.*(So2BF./(Ko2aob+So2BF)).*faob - baob*naob.*((Ko2aob./(Ko2aob+So2BF)).*(Sno2BF+Sno3BF)./(Kno3aob+Sno2BF+Sno3BF)).*faob;
munb=u2.*(Sno2BF./(Kno2nb+Sno2BF)).*(So2BF/(Ko2nb+So2BF)).*fnb - bnb.*(So2BF./(Ko2nb+So2BF)).*fnb - bnb*nnb.*((Ko2nb./(Ko2nb+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nb+Sno2BF+Sno3BF)).*fnb;
munsp=u3.*(Sno2BF./(Kno2nsp+Sno2BF)).*(So2BF/(Ko2nsp+So2BF)).*fnsp - bnsp.*(So2BF./(Ko2nsp+So2BF)).*fnsp - bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nsp+Sno2BF+Sno3BF)).*fnsp;
muamx=u5.*((Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF))).*famx - bamx.*(So2BF./(Ko2amx+So2BF)).*famx - bamx*namx.*((Ko2amx/(Ko2amx+So2BF)).*(Sno2BF+Sno3BF)./(Kno3amx+Sno2BF+Sno3BF)).*famx;
muhan=(u6.*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)) + u6.*(SseBF/(KSsehan+SseBF)).*(Sno3BF/(Kno3han+Sno3BF)).*(Ko2han/(Ko2han+So2BF))+u6.*((SseBF/(KSsehan+SseBF)).*(So2BF/(Ko2han+So2BF))).*fhan)- bhan.*(So2BF/(Ko2han+So2BF)).*fhan - bhan*nhan.*((Ko2han./(Ko2han+So2BF)).*(Sno2BF+Sno3BF)./(Kno3han+Sno2BF+Sno3BF)).*fhan;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L.*muaver+sigma;
VL=VR-AL.*(L + LL); 
%BC-1
SseBFmin=SseBF(nx-(nx-1)); Sno3BFmin=Sno3BF(nx-(nx-1)); Sno2BFmin=Sno2BF(nx-(nx-1)); Snh4BFmin=Snh4BF(nx-(nx-1)); 
So2BFmin=So2BF(nx-(nx-1)); 
%BC-2
SseBFmax=SseBF(nx-1)+L.*DLSse.*(SseBL(nx)-SseBF(nx))./LL; Sno3BFmax=Sno3BF(nx-1)+L.*DLSno3.*(Sno3BL(nx)-Sno3BF(nx))./LL; 
Sno2BFmax=Sno2BF(nx-1)+L.*DLSno2.*(Sno2BL(nx)-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx-1)+L.*DLSnh4.*(Snh4BL(nx)-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+L.*DLSo2.*(So2BL(nx)-So2BF(nx))./LL;  
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
%dSseBLdt=zeros(nx,1);dSno3BLdt=zeros(nx,1);dSno2BLdt=zeros(nx,1);dSnh4BLdt=zeros(nx,1);dSo2BLdt=zeros(nx,1);dXsdt=zeros(nx,1);
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%dLdt=zeros(nx,1);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx-1
    hz=nx-(nx-1);
    z=0.012; zeta=z./L;
    U=L.*muaver*zeta;
    dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
    dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx));
    dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
    dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
    dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
    dXsdt=Qo.*(Xso-Xs)./VL - KH.*((Xs./fhan(nx))./(KX+(Xs/fhan(nx)))).*fhan(nx);
    dSseBFdt(ix)=DSse.*(SseBF(ix+1)-2.*SseBF(ix)+(SseBF(ix-1)))./(hz^2.*L(ix)^2) +zeta.*UL.*(SseBF(ix+1)-SseBF(ix))./L(ix) - (u6*nhan.*(1/YNo2han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) -(u6*nhan.*(1/YNo3han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - (u6/Yhaer).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSno3BFdt(ix)=DSno3.*(Sno3BF(ix+1)-2.*Sno3BF(ix)+(Sno3BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Sno3BF(ix+1)-Sno3BF(ix))./L - (u6*nhan*((1-YNo3han)/(2.86*YNo3han)).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) + (u5/1.14).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix))).*famx(ix) + (u2/Ynb*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) + ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix)) - ((1-fI)/2.86)*baob*naob.*(Ko2aob./(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix)).*faob(ix) - ((1-fI)/2.86)*bnb*nnb.*(Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nb+Sno2BF(ix)+Sno3BF(ix)).*fnb(ix) - ((1-fI)/2.86)*bnsp*nnsp.*(Ko2nsp/(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix)).*fnsp(ix) -  ((1-fI)/2.86)*bamx*namx.*(Ko2amx/(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix)).*famx(ix)-((1-fI)/2.86)*bhan*nhan.*(Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix)).*fhan(ix);
    dSno2BFdt(ix)=DSno2.*(Sno2BF(ix+1)-2.*Sno2BF(ix)+(Sno2BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Sno2BF(ix+1)-Sno2BF(ix))./L - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix))).*fhan(ix)) - (((u5/Yamx)+(u5/1.14)).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)/(Kno2amx+Sno2BF(ix))).*(Ko2amx/(Ko2amx+So2BF(ix)))).*famx(ix) + ((u1/Yaob).*(So2BF(ix)/(Ko2aob+So2(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - ((u2/Ynb).*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) - ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix);
    dSnh4BFdt(ix)=DSnh4.*(Snh4BF(ix+1)-2.*Snh4BF(ix)+(Snh4BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Snh4BF(ix+1)-Snh4BF(ix))./L - (u1*(INBM+1/Yaob).*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) -(u5*(INBM+1/Yamx).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix)))).*famx(ix) - (u2*INBM.*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*(ix) -(u3*INBM.*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - INBM*u6*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) + (INBM-(fI*INXI))*baob*naob.*((Ko2aob/(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb*nnb.*((Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx*namx.*((Ko2amx./(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan*nhan.*((Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix))).*fhan(ix) +(INBM-(fI*INXI))*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSo2BFdt(ix)=DSo2.*(So2BF(ix+1)-2.*So2BF(ix)+(So2BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(So2BF(ix+1)-So2BF(ix))./L - -(((1-Yhaer)/Yhaer)*u6.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) -(((3.43-Yaob)/Yaob)*u1.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - (((1.14-Ynb)/Ynb)*u2.*(Sno2BF(ix)/(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*fnb(ix) - (((1.14-Ynsp)/Ynsp)*u3.*(Sno2BF(ix)/(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - (1-fI)*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) - (1-fI)*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb -(1-fI)*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) - (1-fI)*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) - (1-fI)*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dfaobdt(ix)=(muaob-muaver).*faob(ix) - (U(ix)-zeta.*UL).*(faob(ix+1)-faob(ix))./hz;
    dfnbdt(ix)=(munb-muaver).*fnb(ix) - (U(ix)-zeta.*UL).*(fnb(ix+1)-fnb(ix))./hz;
    dfnspdt(ix)=(munsp-muaver).*fnsp(ix) - (U(ix)-zeta.*UL).*(fnsp(ix+1)-fnsp(ix))./hz;
    dfamxdt(ix)=(muamx-muaver).*famx(ix) - (U(ix)-zeta.*UL).*(famx(ix+1)-famx(ix))./hz;
    dfhandt(ix)=(muhan-muaver).*fhan(ix) - (U(ix)-zeta.*UL).*(fhan(ix+1)-fahan(ix))./hz;
    dLdt=L(ix).*muaver.*zeta+sigma; 
end
    function fval=BC(t,yBC,muaob,munb,munsp,muamx,muhan,muaver)
        fval=zeros(length(yBC),1);
        faobBC=yBC(1,:);
        fnbBC=yBC(2,:);
        fnspBC=yBC(3,:);
        famxBC=yBC(4,:);
        fhanBC=yBC(5,:);
        fval=[(muaob-muaver)*faobBC
            (munb-muaver)*fnbBC
            (munsp-muaver)*fnspBC
            (muamx-muaver)*famxBC
            (muhan-muaver)*fhanBC];
        fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
    end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end

[Torsten]

I cannot run your code since all constants and the call to ode15s is missing.

[KIPROTICH KOSGEY]

sorry for late response but my computer had crashed. Below is the code:

tstart = 0;  
tend = 180*24;   
nt = 14400000;   
xstart = 0;  
xend = 0.012;  
nx = 2000;
t = linspace(tstart,tend,nt);  
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y01 = b(1,:).';  y02=b(2,:).'; y03=b(3,:).';y04=b(4,:).'; y05 = b(5,:).'; y06 = b(6,:).';
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = b(17,:).';

y1=[y01; y02;y03;y04; y05;y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];   
%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
Qo=10; Snh4o=5; KaL=80;
%% execution 
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'InitialStep',1e-10, 'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL),t,y1, options);

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

Error in solution>fun (line 100)
dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));

Error in solution>@(t,y)fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL) (line 48)
[tSol, ySol]=ode15s(@(t,y) fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL),t,y1, 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);

%% feed parameters
function dydt=fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL)
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
dydt=@fval;
muaob=u1.*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF)).*faob - baob.*(So2BF./(Ko2aob+So2BF)).*faob - baob*naob.*((Ko2aob./(Ko2aob+So2BF)).*(Sno2BF+Sno3BF)./(Kno3aob+Sno2BF+Sno3BF)).*faob;
munb=u2.*(Sno2BF./(Kno2nb+Sno2BF)).*(So2BF/(Ko2nb+So2BF)).*fnb - bnb.*(So2BF./(Ko2nb+So2BF)).*fnb - bnb*nnb.*((Ko2nb./(Ko2nb+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nb+Sno2BF+Sno3BF)).*fnb;
munsp=u3.*(Sno2BF./(Kno2nsp+Sno2BF)).*(So2BF/(Ko2nsp+So2BF)).*fnsp - bnsp.*(So2BF./(Ko2nsp+So2BF)).*fnsp - bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nsp+Sno2BF+Sno3BF)).*fnsp;
muamx=u5.*((Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF))).*famx - bamx.*(So2BF./(Ko2amx+So2BF)).*famx - bamx*namx.*((Ko2amx/(Ko2amx+So2BF)).*(Sno2BF+Sno3BF)./(Kno3amx+Sno2BF+Sno3BF)).*famx;
muhan=(u6.*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)) + u6.*(SseBF/(KSsehan+SseBF)).*(Sno3BF/(Kno3han+Sno3BF)).*(Ko2han/(Ko2han+So2BF))+u6.*((SseBF/(KSsehan+SseBF)).*(So2BF/(Ko2han+So2BF))).*fhan)- bhan.*(So2BF/(Ko2han+So2BF)).*fhan - bhan*nhan.*((Ko2han./(Ko2han+So2BF)).*(Sno2BF+Sno3BF)./(Kno3han+Sno2BF+Sno3BF)).*fhan;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L.*muaver+sigma;
VL=VR-AL.*(L + LL); 
%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); 
So2BFmin=So2BF(1); 
%BC-2
SseBFmax=SseBF(nx)+2*(nx-(nx-1))*L.*DLSse.*(SseBL-SseBF(nx))./LL; 
Sno3BFmax=Sno3BF(nx)+2*(nx-(nx-1))*L.*DLSno3.*(Sno3BL-Sno3BF(nx))./LL; 
Sno2BFmax=Sno2BF(nx)+2*(nx-(nx-1))*L.*DLSno2.*(Sno2BL-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx)+2*(nx-(nx-1))*L.*DLSnh4.*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx)+2*(nx-(nx-1))*L.*DLSo2.*(So2BL-So2BF(nx))./LL;  
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
%dSseBLdt=zeros(nx,1);dSno3BLdt=zeros(nx,1);dSno2BLdt=zeros(nx,1);dSnh4BLdt=zeros(nx,1);dSo2BLdt=zeros(nx,1);dXsdt=zeros(nx,1);
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%dLdt=zeros(nx,1);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx
    hz=nx-(nx-1);
    z=0.012; zeta=z./L;
    U=L.*muaver*zeta;
    dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
    dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx));
    dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
    dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
    dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
    dXsdt=Qo.*(Xso-Xs)./VL - KH.*((Xs./fhan(nx))./(KX+(Xs/fhan(nx)))).*fhan(nx);
    dSseBFdt(ix)=DSse.*(SseBF(ix+1)-2.*SseBF(ix)+(SseBF(ix-1)))./(hz^2.*L(ix)^2) +zeta(ix).*UL(ix).*(SseBF(ix+1)-SseBF(ix))./L(ix) - (u6*nhan.*(1/YNo2han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) -(u6*nhan.*(1/YNo3han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - (u6/Yhaer).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSno3BFdt(ix)=DSno3.*(Sno3BF(ix+1)-2.*Sno3BF(ix)+(Sno3BF(ix-1)))./(hz^2.*L(ix)^2) +zeta(ix).*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))./L(ix) - (u6*nhan*((1-YNo3han)/(2.86*YNo3han)).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) + (u5/1.14).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix))).*famx(ix) + (u2/Ynb*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) + ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix)) - ((1-fI)/2.86)*baob*naob.*(Ko2aob./(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix)).*faob(ix) - ((1-fI)/2.86)*bnb*nnb.*(Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nb+Sno2BF(ix)+Sno3BF(ix)).*fnb(ix) - ((1-fI)/2.86)*bnsp*nnsp.*(Ko2nsp/(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix)).*fnsp(ix) -  ((1-fI)/2.86)*bamx*namx.*(Ko2amx/(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix)).*famx(ix)-((1-fI)/2.86)*bhan*nhan.*(Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix)).*fhan(ix);
    dSno2BFdt(ix)=DSno2.*(Sno2BF(ix+1)-2.*Sno2BF(ix)+(Sno2BF(ix-1)))./(hz^2.*L(ix)^2) +zeta(ix).*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))./L(ix) - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix))).*fhan(ix)) - (((u5/Yamx)+(u5/1.14)).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)/(Kno2amx+Sno2BF(ix))).*(Ko2amx/(Ko2amx+So2BF(ix)))).*famx(ix) + ((u1/Yaob).*(So2BF(ix)/(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - ((u2/Ynb).*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) - ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix);
    dSnh4BFdt(ix)=DSnh4.*(Snh4BF(ix+1)-2.*Snh4BF(ix)+(Snh4BF(ix-1)))./(hz^2.*L(ix)^2) +zeta(ix).*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))./L(ix) - (u1*(INBM+1/Yaob).*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) -(u5*(INBM+1/Yamx).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix)))).*famx(ix) - (u2*INBM.*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*(ix) -(u3*INBM.*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - INBM*u6*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) + (INBM-(fI*INXI))*baob*naob.*((Ko2aob/(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb*nnb.*((Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx*namx.*((Ko2amx./(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan*nhan.*((Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix))).*fhan(ix) +(INBM-(fI*INXI))*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSo2BFdt(ix)=DSo2.*(So2BF(ix+1)-2.*So2BF(ix)+(So2BF(ix-1)))./(hz^2.*L(ix)^2) +zeta(ix).*UL(ix).*(So2BF(ix+1)-So2BF(ix))./L(ix) -(((1-Yhaer)/Yhaer)*u6.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) -(((3.43-Yaob)/Yaob)*u1.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - (((1.14-Ynb)/Ynb)*u2.*(Sno2BF(ix)/(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*fnb(ix) - (((1.14-Ynsp)/Ynsp)*u3.*(Sno2BF(ix)/(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - (1-fI)*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) - (1-fI)*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) -(1-fI)*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) - (1-fI)*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) - (1-fI)*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta(ix).*UL(ix)).*(faob(ix+1)-faob(ix))./hz;
    dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta(ix).*UL(ix)).*(fnb(ix+1)-fnb(ix))./hz;
    dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta(ix).*UL(ix)).*(fnsp(ix+1)-fnsp(ix))./hz;
    dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta(ix).*UL(ix)).*(famx(ix+1)-famx(ix))./hz;
    dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta(ix).*UL(ix)).*(fhan(ix+1)-fhan(ix))./hz;
    dLdt=L(ix).*muaver(ix).*zeta(ix)+sigma; 
end
    function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
        fval=zeros(length(yBC),1);
        faobBC=yBC(1:nx,:);
        fnbBC=yBC(nx+1:2*nx,:);
        fnspBC=yBC(2*nx+1:3*nx,:);
        famxBC=yBC(3*nx+1:4*nx,:);
        fhanBC=yBC(4*nx+1:5*nx,:);
        for i=1:nx
            fval=[(muaob-muaver)*faobBC
                (munb-muaver)*fnbBC
                (munsp-muaver)*fnspBC
                (muamx-muaver)*famxBC
                (muhan-muaver)*fhanBC];
        end
        fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
    end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt;faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end

[KIPROTICH KOSGEY]

Dear Torsten,

could you please advise what the problem is in the above code, it is still throwing the index error.

[Walter Roberson]

SseBF=y(7:nx+6,:);

You create SSeBF as a vector (y only has one column)

SseBF=[SseBFmin;SseBFmax];

and there you overwrite the column vector with a different column vector of length 2.

[KIPROTICH KOSGEY]

Thank you Walter!

I have corrected this as :

SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];

However, further down I am getting the error:

Error using vertcat

Dimensions of arrays being concatenated are not consistent.

Error in MBBR>fun (line 155)

Please advise.

[KIPROTICH KOSGEY]

Dear Walter and Torsten,

May I thank you for pointing out errors in my below set of codes. Please I have one last error that I can't figure ot why its coming up yet everything seems fine in my coding: Error using vertcat Dimensions of arrays being concatenated are not consistent. Error in MBBR>fun (line 154)

dydt=[dSseBLdt;dSno3BLdt;......

Edited program:

tstart = 0;  
tend = 180*24;   
nt = 14400000;   
xstart = 0;  
xend = 0.012;  
nx = 2000;
t = linspace(tstart,tend,nt);  
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,6);
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = 1*10^-6;

y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];   
%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
Qo=10; Snh4o=5; KaL=80;
%% execution 
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'InitialStep',1e-10, 'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL),t,y1, options);

Error using vertcat
Dimensions of arrays being concatenated are not consistent.

Error in solution>fun (line 134)
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];

Error in solution>@(t,y)fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL) (line 48)
[tSol, ySol]=ode15s(@(t,y) fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL),t,y1, 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);

%% feed parameters
function dydt=fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL)
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
dydt=@fval;
muaob=u1.*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF)).*faob - baob.*(So2BF./(Ko2aob+So2BF)).*faob - baob*naob.*((Ko2aob./(Ko2aob+So2BF)).*(Sno2BF+Sno3BF)./(Kno3aob+Sno2BF+Sno3BF)).*faob;
munb=u2.*(Sno2BF./(Kno2nb+Sno2BF)).*(So2BF/(Ko2nb+So2BF)).*fnb - bnb.*(So2BF./(Ko2nb+So2BF)).*fnb - bnb*nnb.*((Ko2nb./(Ko2nb+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nb+Sno2BF+Sno3BF)).*fnb;
munsp=u3.*(Sno2BF./(Kno2nsp+Sno2BF)).*(So2BF/(Ko2nsp+So2BF)).*fnsp - bnsp.*(So2BF./(Ko2nsp+So2BF)).*fnsp - bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nsp+Sno2BF+Sno3BF)).*fnsp;
muamx=u5.*((Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF))).*famx - bamx.*(So2BF./(Ko2amx+So2BF)).*famx - bamx*namx.*((Ko2amx/(Ko2amx+So2BF)).*(Sno2BF+Sno3BF)./(Kno3amx+Sno2BF+Sno3BF)).*famx;
muhan=(u6.*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)) + u6.*(SseBF/(KSsehan+SseBF)).*(Sno3BF/(Kno3han+Sno3BF)).*(Ko2han/(Ko2han+So2BF))+u6.*((SseBF/(KSsehan+SseBF)).*(So2BF/(Ko2han+So2BF))).*fhan)- bhan.*(So2BF/(Ko2han+So2BF)).*fhan - bhan*nhan.*((Ko2han./(Ko2han+So2BF)).*(Sno2BF+Sno3BF)./(Kno3han+Sno2BF+Sno3BF)).*fhan;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L.*muaver+sigma;
VL=VR-AL.*(L + LL); 
%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); 
So2BFmin=So2BF(1); 
%BC-2
SseBFmax=SseBF(nx)+2*(nx-(nx-1))*L.*DLSse.*(SseBL-SseBF(nx))./LL; 
Sno3BFmax=Sno3BF(nx)+2*(nx-(nx-1))*L.*DLSno3.*(Sno3BL-Sno3BF(nx))./LL; 
Sno2BFmax=Sno2BF(nx)+2*(nx-(nx-1))*L.*DLSno2.*(Sno2BL-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx)+2*(nx-(nx-1))*L.*DLSnh4.*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx)+2*(nx-(nx-1))*L.*DLSo2.*(So2BL-So2BF(nx))./LL;  
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
dSseBLdt=zeros(nx,1);dSno3BLdt=zeros(nx,1);dSno2BLdt=zeros(nx,1);dSnh4BLdt=zeros(nx,1);dSo2BLdt=zeros(nx,1);dXsdt=zeros(nx,1);
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%dLdt=zeros(nx,1);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx-1
    hz=nx-(nx-1);
    z=0.012; zeta=z./L;
    U=L.*muaver*zeta;
    dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
    dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx));
    dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
    dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
    dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
    dXsdt=Qo.*(Xso-Xs)./VL - KH.*((Xs./fhan(nx))./(KX+(Xs/fhan(nx)))).*fhan(nx);
    dSseBFdt(ix)=DSse.*(SseBF(ix+1)-2.*SseBF(ix)+(SseBF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(SseBF(ix+1)-SseBF(ix))./L - (u6*nhan.*(1/YNo2han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) -(u6*nhan.*(1/YNo3han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - (u6/Yhaer).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSno3BFdt(ix)=DSno3.*(Sno3BF(ix+1)-2.*Sno3BF(ix)+(Sno3BF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))./L - (u6*nhan*((1-YNo3han)/(2.86*YNo3han)).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) + (u5/1.14).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix))).*famx(ix) + (u2/Ynb*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) + ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix)) - ((1-fI)/2.86)*baob*naob.*(Ko2aob./(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix)).*faob(ix) - ((1-fI)/2.86)*bnb*nnb.*(Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nb+Sno2BF(ix)+Sno3BF(ix)).*fnb(ix) - ((1-fI)/2.86)*bnsp*nnsp.*(Ko2nsp/(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix)).*fnsp(ix) -  ((1-fI)/2.86)*bamx*namx.*(Ko2amx/(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix)).*famx(ix)-((1-fI)/2.86)*bhan*nhan.*(Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix)).*fhan(ix);
    dSno2BFdt(ix)=DSno2.*(Sno2BF(ix+1)-2.*Sno2BF(ix)+(Sno2BF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))./L - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix))).*fhan(ix)) - (((u5/Yamx)+(u5/1.14)).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)/(Kno2amx+Sno2BF(ix))).*(Ko2amx/(Ko2amx+So2BF(ix)))).*famx(ix) + ((u1/Yaob).*(So2BF(ix)/(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - ((u2/Ynb).*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) - ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix);
    dSnh4BFdt(ix)=DSnh4.*(Snh4BF(ix+1)-2.*Snh4BF(ix)+(Snh4BF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))./L - (u1*(INBM+1/Yaob).*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) -(u5*(INBM+1/Yamx).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix)))).*famx(ix) - (u2*INBM.*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*(ix) -(u3*INBM.*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - INBM*u6*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) + (INBM-(fI*INXI))*baob*naob.*((Ko2aob/(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb*nnb.*((Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx*namx.*((Ko2amx./(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan*nhan.*((Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix))).*fhan(ix) +(INBM-(fI*INXI))*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSo2BFdt(ix)=DSo2.*(So2BF(ix+1)-2.*So2BF(ix)+(So2BF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(So2BF(ix+1)-So2BF(ix))./L -(((1-Yhaer)/Yhaer)*u6.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) -(((3.43-Yaob)/Yaob)*u1.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - (((1.14-Ynb)/Ynb)*u2.*(Sno2BF(ix)/(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*fnb(ix) - (((1.14-Ynsp)/Ynsp)*u3.*(Sno2BF(ix)/(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - (1-fI)*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) - (1-fI)*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) -(1-fI)*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) - (1-fI)*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) - (1-fI)*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta.*UL(ix)).*(faob(ix+1)-faob(ix))./hz;
    dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta.*UL(ix)).*(fnb(ix+1)-fnb(ix))./hz;
    dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta.*UL(ix)).*(fnsp(ix+1)-fnsp(ix))./hz;
    dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta.*UL(ix)).*(famx(ix+1)-famx(ix))./hz;
    dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta.*UL(ix)).*(fhan(ix+1)-fhan(ix))./hz;
    dLdt=L.*muaver(ix).*zeta+sigma; 
end
    function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
        fval=zeros(length(yBC),1);
        faobBC=yBC(1,:);
        fnbBC=yBC(2,:);
        fnspBC=yBC(3,:);
        famxBC=yBC(4,:);
        fhanBC=yBC(5,:);
        for i=1
            fval=[(muaob(i)-muaver(i))*faobBC(i)
                (munb(i)-muaver(i))*fnbBC(i)
                (munsp(i)-muaver(i))*fnspBC(i)
                (muamx(i)-muaver(i))*famxBC(i)
                (muhan(i)-muaver(i))*fhanBC(i)];
        end
        fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
    end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end

[Walter Roberson]

Observe that dSseBLdt is 2000 x 2000.

It is likely that you have accidentally combined row vector calculations with column vector calculations. For example, (1:3)' + [10 20 30] will not result in [11 22 33] and will instead give a 3 x 3 result.

tstart = 0;  
tend = 180*24;   
nt = 14400000;   
xstart = 0;  
xend = 0.012;  
nx = 2000;
t = linspace(tstart,tend,nt);  
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,6);
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = 1*10^-6;

y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];   
%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
Qo=10; Snh4o=5; KaL=80;
%% execution 
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'InitialStep',1e-10, 'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL),t,y1, options);

  Name               Size                 Bytes  Class              Attributes

  AL                 1x1                      8  double                       
  DLSnh4             1x1                      8  double                       
  DLSno2             1x1                      8  double                       
  DLSno3             1x1                      8  double                       
  DLSo2              1x1                      8  double                       
  DLSse              1x1                      8  double                       
  DSnh4              1x1                      8  double                       
  DSno2              1x1                      8  double                       
  DSno3              1x1                      8  double                       
  DSo2               1x1                      8  double                       
  DSse               1x1                      8  double                       
  INBM               1x1                      8  double                       
  INXI               1x1                      8  double                       
  KH                 1x1                      8  double                       
  KSsehan            1x1                      8  double                       
  KX                 1x1                      8  double                       
  KaL                1x1                      8  double                       
  Knh4amx            1x1                      8  double                       
  Knh4aob            1x1                      8  double                       
  Kno2amx            1x1                      8  double                       
  Kno2han            1x1                      8  double                       
  Kno2nb             1x1                      8  double                       
  Kno2nsp            1x1                      8  double                       
  Kno3amx            1x1                      8  double                       
  Kno3aob            1x1                      8  double                       
  Kno3han            1x1                      8  double                       
  Kno3nb             1x1                      8  double                       
  Kno3nsp            1x1                      8  double                       
  Ko2amx             1x1                      8  double                       
  Ko2aob             1x1                      8  double                       
  Ko2han             1x1                      8  double                       
  Ko2nb              1x1                      8  double                       
  Ko2nsp             1x1                      8  double                       
  L                  1x1                      8  double                       
  LL                 1x1                      8  double                       
  Qo                 1x1                      8  double                       
  Snh4BF          2001x1                  16008  double                       
  Snh4BFmax          1x1                      8  double                       
  Snh4BFmin          1x1                      8  double                       
  Snh4BL             1x1                      8  double                       
  Snh4o              1x1                      8  double                       
  Sno2BF          2001x1                  16008  double                       
  Sno2BFmax          1x1                      8  double                       
  Sno2BFmin          1x1                      8  double                       
  Sno2BL             1x1                      8  double                       
  Sno2o              1x1                      8  double                       
  Sno3BF          2001x1                  16008  double                       
  Sno3BFmax          1x1                      8  double                       
  Sno3BFmin          1x1                      8  double                       
  Sno3BL             1x1                      8  double                       
  Sno3o              1x1                      8  double                       
  So2BF           2001x1                  16008  double                       
  So2BFmax           1x1                      8  double                       
  So2BFmin           1x1                      8  double                       
  So2BL              1x1                      8  double                       
  So2o               1x1                      8  double                       
  SseBF           2001x1                  16008  double                       
  SseBFmax           1x1                      8  double                       
  SseBFmin           1x1                      8  double                       
  SseBL              1x1                      8  double                       
  Sseo               1x1                      8  double                       
  U               2000x2000            32000000  double                       
  UL              2000x2000            32000000  double                       
  VL                 1x1                      8  double                       
  VR                 1x1                      8  double                       
  Xs                 1x1                      8  double                       
  Xso                1x1                      8  double                       
  YNo2han            1x1                      8  double                       
  YNo3han            1x1                      8  double                       
  Yamx               1x1                      8  double                       
  Yaob               1x1                      8  double                       
  Yhaer              1x1                      8  double                       
  Ynb                1x1                      8  double                       
  Ynsp               1x1                      8  double                       
  bamx               1x1                      8  double                       
  baob               1x1                      8  double                       
  bhan               1x1                      8  double                       
  bnb                1x1                      8  double                       
  bnsp               1x1                      8  double                       
  dLdt               1x1                      8  double                       
  dSnh4BFdt       2000x1                  16000  double                       
  dSnh4BLdt       2000x2000            32000000  double                       
  dSno2BFdt       2000x1                  16000  double                       
  dSno2BLdt       2000x2000            32000000  double                       
  dSno3BFdt       2000x1                  16000  double                       
  dSno3BLdt       2000x2000            32000000  double                       
  dSo2BFdt        2000x1                  16000  double                       
  dSo2BLdt        2000x2000            32000000  double                       
  dSseBFdt        2000x1                  16000  double                       
  dSseBLdt        2000x2000            32000000  double                       
  dXsdt              1x1                      8  double                       
  dfamxdt         2000x1                  16000  double                       
  dfaobdt         2000x1                  16000  double                       
  dfhandt         2000x1                  16000  double                       
  dfnbdt          2000x1                  16000  double                       
  dfnspdt         2000x1                  16000  double                       
  dydt               1x1                     32  function_handle              
  fI                 1x1                      8  double                       
  famx            2001x1                  16008  double                       
  famxBC             1x1                      8  double                       
  faob            2001x1                  16008  double                       
  faobBC             1x1                      8  double                       
  fhan            2001x1                  16008  double                       
  fhanBC             1x1                      8  double                       
  fnb             2001x1                  16008  double                       
  fnbBC              1x1                      8  double                       
  fnsp            2001x1                  16008  double                       
  fnspBC             1x1                      8  double                       
  hz                 1x1                      8  double                       
  ix                 1x1                      8  double                       
  muamx           2000x2000            32000000  double                       
  muaob           2000x1                  16000  double                       
  muaver          2000x2000            32000000  double                       
  muhan           2000x2000            32000000  double                       
  munb            2000x2000            32000000  double                       
  munsp           2000x2000            32000000  double                       
  namx               1x1                      8  double                       
  naob               1x1                      8  double                       
  nhan               1x1                      8  double                       
  nnb                1x1                      8  double                       
  nnsp               1x1                      8  double                       
  nx                 1x1                      8  double                       
  sigma              1x1                      8  double                       
  t                  1x1                      8  double                       
  u1                 1x1                      8  double                       
  u2                 1x1                      8  double                       
  u3                 1x1                      8  double                       
  u5                 1x1                      8  double                       
  u6                 1x1                      8  double                       
  y              20007x1                 160056  double                       
  z                  1x1                      8  double                       
  zeta               1x1                      8  double                       

Error using vertcat
Dimensions of arrays being concatenated are not consistent.

Error in solution>fun (line 135)
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];

Error in solution>@(t,y)fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL) (line 48)
[tSol, ySol]=ode15s(@(t,y) fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL),t,y1, 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);

%% feed parameters
function dydt=fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL)
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
dydt=@fval;
muaob=u1.*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF)).*faob - baob.*(So2BF./(Ko2aob+So2BF)).*faob - baob*naob.*((Ko2aob./(Ko2aob+So2BF)).*(Sno2BF+Sno3BF)./(Kno3aob+Sno2BF+Sno3BF)).*faob;
munb=u2.*(Sno2BF./(Kno2nb+Sno2BF)).*(So2BF/(Ko2nb+So2BF)).*fnb - bnb.*(So2BF./(Ko2nb+So2BF)).*fnb - bnb*nnb.*((Ko2nb./(Ko2nb+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nb+Sno2BF+Sno3BF)).*fnb;
munsp=u3.*(Sno2BF./(Kno2nsp+Sno2BF)).*(So2BF/(Ko2nsp+So2BF)).*fnsp - bnsp.*(So2BF./(Ko2nsp+So2BF)).*fnsp - bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nsp+Sno2BF+Sno3BF)).*fnsp;
muamx=u5.*((Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF))).*famx - bamx.*(So2BF./(Ko2amx+So2BF)).*famx - bamx*namx.*((Ko2amx/(Ko2amx+So2BF)).*(Sno2BF+Sno3BF)./(Kno3amx+Sno2BF+Sno3BF)).*famx;
muhan=(u6.*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)) + u6.*(SseBF/(KSsehan+SseBF)).*(Sno3BF/(Kno3han+Sno3BF)).*(Ko2han/(Ko2han+So2BF))+u6.*((SseBF/(KSsehan+SseBF)).*(So2BF/(Ko2han+So2BF))).*fhan)- bhan.*(So2BF/(Ko2han+So2BF)).*fhan - bhan*nhan.*((Ko2han./(Ko2han+So2BF)).*(Sno2BF+Sno3BF)./(Kno3han+Sno2BF+Sno3BF)).*fhan;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L.*muaver+sigma;
VL=VR-AL.*(L + LL); 
%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); 
So2BFmin=So2BF(1); 
%BC-2
SseBFmax=SseBF(nx)+2*(nx-(nx-1))*L.*DLSse.*(SseBL-SseBF(nx))./LL; 
Sno3BFmax=Sno3BF(nx)+2*(nx-(nx-1))*L.*DLSno3.*(Sno3BL-Sno3BF(nx))./LL; 
Sno2BFmax=Sno2BF(nx)+2*(nx-(nx-1))*L.*DLSno2.*(Sno2BL-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx)+2*(nx-(nx-1))*L.*DLSnh4.*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx)+2*(nx-(nx-1))*L.*DLSo2.*(So2BL-So2BF(nx))./LL;  
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
dSseBLdt=zeros(nx,1);dSno3BLdt=zeros(nx,1);dSno2BLdt=zeros(nx,1);dSnh4BLdt=zeros(nx,1);dSo2BLdt=zeros(nx,1);dXsdt=zeros(nx,1);
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%dLdt=zeros(nx,1);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx-1
    hz=nx-(nx-1);
    z=0.012; zeta=z./L;
    U=L.*muaver*zeta;
    dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
    dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx));
    dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
    dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
    dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
    dXsdt=Qo.*(Xso-Xs)./VL - KH.*((Xs./fhan(nx))./(KX+(Xs/fhan(nx)))).*fhan(nx);
    dSseBFdt(ix)=DSse.*(SseBF(ix+1)-2.*SseBF(ix)+(SseBF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(SseBF(ix+1)-SseBF(ix))./L - (u6*nhan.*(1/YNo2han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) -(u6*nhan.*(1/YNo3han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - (u6/Yhaer).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSno3BFdt(ix)=DSno3.*(Sno3BF(ix+1)-2.*Sno3BF(ix)+(Sno3BF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))./L - (u6*nhan*((1-YNo3han)/(2.86*YNo3han)).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) + (u5/1.14).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix))).*famx(ix) + (u2/Ynb*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) + ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix)) - ((1-fI)/2.86)*baob*naob.*(Ko2aob./(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix)).*faob(ix) - ((1-fI)/2.86)*bnb*nnb.*(Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nb+Sno2BF(ix)+Sno3BF(ix)).*fnb(ix) - ((1-fI)/2.86)*bnsp*nnsp.*(Ko2nsp/(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix)).*fnsp(ix) -  ((1-fI)/2.86)*bamx*namx.*(Ko2amx/(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix)).*famx(ix)-((1-fI)/2.86)*bhan*nhan.*(Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix)).*fhan(ix);
    dSno2BFdt(ix)=DSno2.*(Sno2BF(ix+1)-2.*Sno2BF(ix)+(Sno2BF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))./L - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix))).*fhan(ix)) - (((u5/Yamx)+(u5/1.14)).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)/(Kno2amx+Sno2BF(ix))).*(Ko2amx/(Ko2amx+So2BF(ix)))).*famx(ix) + ((u1/Yaob).*(So2BF(ix)/(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - ((u2/Ynb).*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) - ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix);
    dSnh4BFdt(ix)=DSnh4.*(Snh4BF(ix+1)-2.*Snh4BF(ix)+(Snh4BF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))./L - (u1*(INBM+1/Yaob).*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) -(u5*(INBM+1/Yamx).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix)))).*famx(ix) - (u2*INBM.*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*(ix) -(u3*INBM.*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - INBM*u6*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) + (INBM-(fI*INXI))*baob*naob.*((Ko2aob/(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb*nnb.*((Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx*namx.*((Ko2amx./(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan*nhan.*((Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix))).*fhan(ix) +(INBM-(fI*INXI))*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dSo2BFdt(ix)=DSo2.*(So2BF(ix+1)-2.*So2BF(ix)+(So2BF(ix-1)))./(hz^2.*L^2) +zeta.*UL(ix).*(So2BF(ix+1)-So2BF(ix))./L -(((1-Yhaer)/Yhaer)*u6.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) -(((3.43-Yaob)/Yaob)*u1.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - (((1.14-Ynb)/Ynb)*u2.*(Sno2BF(ix)/(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*fnb(ix) - (((1.14-Ynsp)/Ynsp)*u3.*(Sno2BF(ix)/(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - (1-fI)*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) - (1-fI)*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) -(1-fI)*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) - (1-fI)*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) - (1-fI)*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
    dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta.*UL(ix)).*(faob(ix+1)-faob(ix))./hz;
    dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta.*UL(ix)).*(fnb(ix+1)-fnb(ix))./hz;
    dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta.*UL(ix)).*(fnsp(ix+1)-fnsp(ix))./hz;
    dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta.*UL(ix)).*(famx(ix+1)-famx(ix))./hz;
    dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta.*UL(ix)).*(fhan(ix+1)-fhan(ix))./hz;
    dLdt=L.*muaver(ix).*zeta+sigma; 
end
    function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
        fval=zeros(length(yBC),1);
        faobBC=yBC(1,:);
        fnbBC=yBC(2,:);
        fnspBC=yBC(3,:);
        famxBC=yBC(4,:);
        fhanBC=yBC(5,:);
        for i=1
            fval=[(muaob(i)-muaver(i))*faobBC(i)
                (munb(i)-muaver(i))*fnbBC(i)
                (munsp(i)-muaver(i))*fnspBC(i)
                (muamx(i)-muaver(i))*famxBC(i)
                (muhan(i)-muaver(i))*fhanBC(i)];
        end
        fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
    end
whos
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end

[KIPROTICH KOSGEY]

Hi Walter,

How do you get the sizes of the variables before the program is executed?

If I could manage to do this then I could be checking their attributes as I make corrections.

[Torsten]

https://uk.mathworks.com/help/matlab/matlab_prog/debugging-process-and-features.html

[KIPROTICH KOSGEY]

Thank you Torsten and Walter.

I have managed to remove all the bugs. But code execution is very slow and takes up a lot of memory. equations are also very stiff.

To save on memory, I have done program execution in bits (timewise). On stiffness, I am re-running the program from the time it fails using the final output from the previous run as input for the new runs. However, program execution is very slow. Could there be a way to speed the execution?

Please see the codes attached.

tstart = 0;  tstart1=tSol(end);

The end operator must be used within an array index expression.

tend = 10*24;   
nt = 80000;   
xstart = 0;  
xend = 0.012;  
nx = 1000;
t = linspace(tstart,tend,nt);  t1 = linspace(tstart1,tend,nt); 
x = linspace(xstart,xend,nx);

% Method-of-lines solution
x = x.';

% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];

a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,6);
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = 1*10^-6;
   
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17]; 

%% execution 
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1);
y2=[ySol(end,1); ySol(end,2);ySol(end,3);ySol(end,4); ySol(end,5);ySol(end,6);ySol(end,7:nx+6)'; ySol(end,nx+7:2*nx+6)';ySol(end,2*nx+7:3*nx+6)';ySol(end,3*nx+7:4*nx+6)'; ySol(end,4*nx+7:5*nx+6)';ySol(end,5*nx+7:6*nx+6)';ySol(end,6*nx+7:7*nx+6)'; ySol(end,7*nx+7:8*nx+6)';ySol(end,8*nx+7:9*nx+6)';ySol(end,9*nx+7:10*nx+6)'; ySol(end,10*nx+7)];
[tSol1, ySol1]=ode15s(@(t1,y) MBBR(t1,y),t1,y2);


writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions

function u0 =  oscic(x)
  u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end

function dydt=MBBR(t,y)

tstart = 0;  
tend = 10*24;   
nt = 80000;   
xstart = 0;  
xend = 0.012;  
nx = 1000;
t = linspace(tstart,tend,nt);  
x = linspace(xstart,xend,nx);

% Method-of-lines solution
x = x.';

% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];

%% constants
T=301;
%aob

O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; n=0.5; 
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;

AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
KaL=0;

%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);


dydt=@fval;

muaob=u1*(So2BF/(Ko2aob+So2BF))*(Snh4BF./(Knh4aob+Snh4BF));                             MDaobaer =baob*(So2BF./(Ko2aob+So2BF));   MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)));
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF));                                MDnbaer= bnb*(So2BF./(Ko2nb+So2BF));      MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF));
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF));                             MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF));  MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF));
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF));  MDamxaer= bamx*(So2BF./(Ko2amx+So2BF));  MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF));
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)); muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF)); muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF)); MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF));
muhan=muhanSo2+muhanSno3+muhanSno2;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;

UL=L*muaver+sigma;
VL=VR-AL*(L + LL); 
z=0.012; zeta=z/L;
U=L*muaver*zeta;
m=0.01; %mass/unit_depth (mg/m)
fill_ratio=0.43;
VC=fill_ratio*VR; %volume occupied by carriers
X=5000;

if and (ge(t,0),le(t,5))
   Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
 elseif  and (gt(t,5),le(t,9))     
     Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
 elseif and (gt(t,9), le(t,12))
     Snh4o=909;Qo=0.132*VL*10^3/(909*24);
 elseif and (gt(t,12), le(t,15))
     Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
 elseif and (gt(t,15), le(t,21))
     Snh4o=937;Qo=0.185*VL*10^3/(937*24);
 elseif and (gt(t,21),le(t,33))
     Snh4o=937;Qo=0.2286*VL*10^3/(937*24); 
 elseif and (gt(t,33),le(t,35))
     Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
 elseif and (gt(t,35),le(t,37))
     Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
 elseif and (gt(t,37),le(t,49))
     Snh4o=937;Qo=0.348*VL*10^3/(937*24);
 elseif and (gt(t,49),le(t,77))
     Snh4o=937;Qo=0.448*VL*10^3/(937*24);
 elseif and (gt(t,77),le(t,91))
     Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
 elseif and (gt(t,91),le(t,119))
     Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
 elseif and (gt(t,119),le(t,127))
     Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
  elseif and (gt(t,127),le(t,134))
     Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
 elseif and (gt(t,134),le(t,140))
     Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
 elseif and (gt(t,140),le(t,148))
     Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
 elseif and (gt(t,148),le(t,155))
     Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
 elseif and (gt(t,155), le(t,161))
     Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
 elseif and (gt(t,161),le(t,169))
     Snh4o=1078;Qo=0.564*VL*10^3/(1078*24);
 else
     Snh4o=928; Qo=0.47*VL*10^3/(928*24);
end

while So2BL<0.5
        KaL=80;
end

%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); 
So2BFmin=So2BF(1); 

%BC-2
SseBFmax=SseBF(nx-1)+2*(nx-(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL; 
Sno3BFmax=Sno3BF(nx-1)+2*(nx-(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL; 
Sno2BFmax=Sno2BF(nx-1)+2*(nx-(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx-1)+2*(nx-(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(nx-(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;  

faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];

dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);



%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx-1
            
    dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
    dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BF(nx));
    dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BF(nx));
    dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BF(nx));
    dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BF(nx));
    dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/fhan(nx))/(KX+(Xs/fhan(nx))))*fhan(nx);
    dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(SseBF(ix+1)-SseBF(ix))/L - n*(1/YNo2han)*muhanSno2(ix)*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix)*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix)*fhan(ix)*X;
    dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Sno3BF(ix+1)-Sno3BF(ix))/L - n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix)*fhan(ix)*X + (1/1.14)*muamx(ix)*famx(ix)*X + (1/Ynb)*munb(ix)*fnb(ix)*X + (1/Ynsp)*munsp(ix)*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix)*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix)*fnb(ix) - ((1-fI)/2.86)*MDnspana(ix)*fnsp(ix)*X -  ((1-fI)/2.86)*MDamxana(ix)*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix)*fhan(ix)*X;
    dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Sno2BF(ix+1)-Sno2BF(ix))/L - ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix)*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix)*famx(ix)*X + (1/Yaob)*muaob(ix)*faob(ix)*X - (1/Ynb)*munb(ix)*fnb(ix)*X - (1/Ynsp)*munb(ix)*fnsp(ix)*X;
    dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Snh4BF(ix+1)-Snh4BF(ix))/L - (INBM+1/Yaob)*muaob(ix)*faob(ix)*X -(INBM+1/Yamx)*muamx(ix)*famx(ix)*X - INBM*munb(ix)*fnb(ix)*X -INBM*munsp(ix)*fnsp(ix)*X - INBM*n*muhanSno2(ix)*fhan(ix)*X - INBM*n*muhanSno3(ix)*fhan(ix)*X - INBM*muhanSo2(ix)*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix)*faob(ix)*X + (INBM-(fI*INXI))*MDnbana(ix)*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix)*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix)*famx(ix)*X + (INBM-(fI*INXI))*MDhanana(ix)*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix)*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix)*fnb(ix)*X +(INBM-(fI*INXI))*MDnspaer(ix)*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix)*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix)*fhan(ix)*X;
    dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L -((1-Yhaer)/Yhaer)*muhanSo2(ix)*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix)*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix)*fnb(ix)*X - ((1.14-Ynsp)/Ynsp)*munsp(ix)*fnsp(ix)*X - (1-fI)*MDaobaer(ix)*faob(ix)*X - (1-fI)*MDnbaer(ix)*fnb(ix)*X -(1-fI)*MDnspaer(ix)*fnsp(ix)*X - (1-fI)*MDamxaer(ix)*famx(ix)*X - (1-fI)*MDhanaer(ix)*fhan(ix);
    dfaobdt(ix)=(muaob(ix)-muaver(ix))*faob(ix) - (U(ix)-zeta*UL(ix))*(faob(ix+1)-faob(ix))/(x(nx+1)-x(nx-1));
    dfnbdt(ix)=(munb(ix)-muaver(ix))*fnb(ix) - (U(ix)-zeta*UL(ix))*(fnb(ix+1)-fnb(ix))/(x(nx+1)-x(nx-1));
    dfnspdt(ix)=(munsp(ix)-muaver(ix))*fnsp(ix) - (U(ix)-zeta*UL(ix))*(fnsp(ix+1)-fnsp(ix))/(x(nx+1)-x(nx-1));
    dfamxdt(ix)=(muamx(ix)-muaver(ix))*famx(ix) - (U(ix)-zeta*UL(ix))*(famx(ix+1)-famx(ix))/(x(nx+1)-x(nx-1));
    dfhandt(ix)=(muhan(ix)-muaver(ix))*fhan(ix) - (U(ix)-zeta*UL(ix))*(fhan(ix+1)-fhan(ix))/(x(nx+1)-x(nx-1));
    dLdt=L*muaver(nx)*zeta+sigma;
        
end
    function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
            fval=zeros(length(yBC),1);
            faobBC=yBC(1,:);
            fnbBC=yBC(2,:);
            fnspBC=yBC(3,:);
            famxBC=yBC(4,:);
            fhanBC=yBC(5,:);
            for i=1
                fval=[(muaob-muaver)*faobBC(i)
                      (munb-muaver)*fnbBC(i)
                      (munsp-muaver)*fnspBC(i)
                      (muamx-muaver)*famxBC(i)
                      (muhan-muaver)*fhanBC(i)];
            end
            fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
        end

dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end

[Torsten]

I can't make the two .m files work.

[KIPROTICH KOSGEY]

this is the code:

tstart = 0;  
tend = 10*24;   
nt = 160000;   
xstart = 0;  
xend = 0.012;  
nx = 500;
t = linspace(tstart,tend,nt);   
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,6);
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = 1*10^-6;

y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17]; 
%% execution 
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1);
%writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 =  oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end

function dydt=MBBR(t,y)
tstart = 0;  
tend = 10*24;   
nt = 160000;   
xstart = 0;  
xend = 0.012;  
nx = 500;
t = linspace(tstart,tend,nt);   
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; n=0.5; 
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
KaL=0;
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
dydt=@fval;
muaob=u1*(So2BF/(Ko2aob+So2BF))*(Snh4BF./(Knh4aob+Snh4BF));                             MDaobaer =baob*(So2BF./(Ko2aob+So2BF));   MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)));
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF));                                MDnbaer= bnb*(So2BF./(Ko2nb+So2BF));      MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF));
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF));                             MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF));  MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF));
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF));  MDamxaer= bamx*(So2BF./(Ko2amx+So2BF));  MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF));
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)); muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF)); muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF)); MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF));
muhan=muhanSo2+muhanSno3+muhanSno2;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L*muaver+sigma;
VL=VR-AL*(L + LL); 
z=0.012; zeta=z/L;
U=L*muaver*zeta;
m=0.01; %mass/unit_depth (mg/m)
fill_ratio=0.43;
VC=fill_ratio*VR; %volume occupied by carriers
X=5000;
if and (ge(t,0),le(t,5))
    Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
elseif  and (gt(t,5),le(t,9))     
    Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
elseif and (gt(t,9), le(t,12))
    Snh4o=909;Qo=0.132*VL*10^3/(909*24);
elseif and (gt(t,12), le(t,15))
    Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
elseif and (gt(t,15), le(t,21))
    Snh4o=937;Qo=0.185*VL*10^3/(937*24);
elseif and (gt(t,21),le(t,33))
    Snh4o=937;Qo=0.2286*VL*10^3/(937*24); 
elseif and (gt(t,33),le(t,35))
    Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
elseif and (gt(t,35),le(t,37))
    Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
elseif and (gt(t,37),le(t,49))
    Snh4o=937;Qo=0.348*VL*10^3/(937*24);
elseif and (gt(t,49),le(t,77))
    Snh4o=937;Qo=0.448*VL*10^3/(937*24);
elseif and (gt(t,77),le(t,91))
    Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
elseif and (gt(t,91),le(t,119))
    Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
elseif and (gt(t,119),le(t,127))
    Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
elseif and (gt(t,127),le(t,134))
    Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
elseif and (gt(t,134),le(t,140))
    Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
elseif and (gt(t,140),le(t,148))
    Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
elseif and (gt(t,148),le(t,155))
    Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
elseif and (gt(t,155), le(t,161))
    Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
elseif and (gt(t,161),le(t,169))
    Snh4o=1078;Qo=0.564*VL*10^3/(1078*24);
else
    Snh4o=928; Qo=0.47*VL*10^3/(928*24);
end
while So2BL<0.5
    KaL=80;
end
%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); 
So2BFmin=So2BF(1); 
%BC-2
SseBFmax=SseBF(nx-1)+2*(nx-(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL; 
Sno3BFmax=Sno3BF(nx-1)+2*(nx-(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL; 
Sno2BFmax=Sno2BF(nx-1)+2*(nx-(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx-1)+2*(nx-(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(nx-(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;  
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx-1
    
    dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
    dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BF(nx));
    dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BF(nx));
    dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BF(nx));
    dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BF(nx));
    dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/fhan(nx))/(KX+(Xs/fhan(nx))))*fhan(nx);
    dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(SseBF(ix+1)-SseBF(ix))/L - n*(1/YNo2han)*muhanSno2(ix)*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix)*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix)*fhan(ix)*X;
    dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Sno3BF(ix+1)-Sno3BF(ix))/L - n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix)*fhan(ix)*X + (1/1.14)*muamx(ix)*famx(ix)*X + (1/Ynb)*munb(ix)*fnb(ix)*X + (1/Ynsp)*munsp(ix)*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix)*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix)*fnb(ix) - ((1-fI)/2.86)*MDnspana(ix)*fnsp(ix)*X -  ((1-fI)/2.86)*MDamxana(ix)*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix)*fhan(ix)*X;
    dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Sno2BF(ix+1)-Sno2BF(ix))/L - ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix)*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix)*famx(ix)*X + (1/Yaob)*muaob(ix)*faob(ix)*X - (1/Ynb)*munb(ix)*fnb(ix)*X - (1/Ynsp)*munb(ix)*fnsp(ix)*X;
    dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Snh4BF(ix+1)-Snh4BF(ix))/L - (INBM+1/Yaob)*muaob(ix)*faob(ix)*X -(INBM+1/Yamx)*muamx(ix)*famx(ix)*X - INBM*munb(ix)*fnb(ix)*X -INBM*munsp(ix)*fnsp(ix)*X - INBM*n*muhanSno2(ix)*fhan(ix)*X - INBM*n*muhanSno3(ix)*fhan(ix)*X - INBM*muhanSo2(ix)*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix)*faob(ix)*X + (INBM-(fI*INXI))*MDnbana(ix)*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix)*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix)*famx(ix)*X + (INBM-(fI*INXI))*MDhanana(ix)*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix)*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix)*fnb(ix)*X +(INBM-(fI*INXI))*MDnspaer(ix)*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix)*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix)*fhan(ix)*X;
    dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L -((1-Yhaer)/Yhaer)*muhanSo2(ix)*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix)*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix)*fnb(ix)*X - ((1.14-Ynsp)/Ynsp)*munsp(ix)*fnsp(ix)*X - (1-fI)*MDaobaer(ix)*faob(ix)*X - (1-fI)*MDnbaer(ix)*fnb(ix)*X -(1-fI)*MDnspaer(ix)*fnsp(ix)*X - (1-fI)*MDamxaer(ix)*famx(ix)*X - (1-fI)*MDhanaer(ix)*fhan(ix);
    dfaobdt(ix)=(muaob(ix)-muaver(ix))*faob(ix) - (U(ix)-zeta*UL(ix))*(faob(ix+1)-faob(ix))/(x(nx+1)-x(nx-1));
    dfnbdt(ix)=(munb(ix)-muaver(ix))*fnb(ix) - (U(ix)-zeta*UL(ix))*(fnb(ix+1)-fnb(ix))/(x(nx+1)-x(nx-1));
    dfnspdt(ix)=(munsp(ix)-muaver(ix))*fnsp(ix) - (U(ix)-zeta*UL(ix))*(fnsp(ix+1)-fnsp(ix))/(x(nx+1)-x(nx-1));
    dfamxdt(ix)=(muamx(ix)-muaver(ix))*famx(ix) - (U(ix)-zeta*UL(ix))*(famx(ix+1)-famx(ix))/(x(nx+1)-x(nx-1));
    dfhandt(ix)=(muhan(ix)-muaver(ix))*fhan(ix) - (U(ix)-zeta*UL(ix))*(fhan(ix+1)-fhan(ix))/(x(nx+1)-x(nx-1));
    dLdt=L*muaver(nx)*zeta+sigma;
    
end
    function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
        fval=zeros(length(yBC),1);
        faobBC=yBC(1,:);
        fnbBC=yBC(2,:);
        fnspBC=yBC(3,:);
        famxBC=yBC(4,:);
        fhanBC=yBC(5,:);
        for i=1
            fval=[(muaob-muaver)*faobBC(i)
                (munb-muaver)*fnbBC(i)
                (munsp-muaver)*fnspBC(i)
                (muamx-muaver)*famxBC(i)
                (muhan-muaver)*fhanBC(i)];
        end
        fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
    end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end

[KIPROTICH KOSGEY]

Dear Torsten and Walter,

Kindly, could there be a way of circumventing the error: 'Warning: Failure at t=1.692955e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

> In ode15s (line 662) In MBBRScript (line 28)' ?

[Torsten]

Did you look at the results obtained so far (for times < 1.7 s) ? Do they make sense ?

Sorry, but it's really impossible to debug your code. Split the long expressions in smaller pieces according to the different parts of your equations, insert comments etc.

Don't treat all your equations in one big loop, but each equation separately.

It makes no sense to compute the ODEs for each ix = 2:nx-1 again and again.

    dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
    dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BF(nx));
    dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BF(nx));
    dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BF(nx));
    dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BF(nx));
    dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/fhan(nx))/(KX+(Xs/fhan(nx))))*fhan(nx); 
    dLdt=L*muaver(nx)*zeta+sigma;

Further I don't see a call to "BC" in the visible part of your code. Maybe it's hidden somewhere in the endless lines where you compute the derivatives.

[KIPROTICH KOSGEY]

Hi Torsten,

I figured out that computing ODEs with SseBF(nx),Sno3BF(nx),Sno2BF(nx),Snh4BF(nx) and So2BF(nx) while ix= 2:nx-1 resulted in a constant values of ODEs since I was using a none-existent values. Therefore I have changed ix=2:nx.

However, there seem to be a problem with the code, as matlab seem to be getting confused since it never finishes computations even after reducing nx and nt to 10 each. I have left MATLAB to run for over 24 hours but still it doesn't finish computing.

I have relooked at the codes and realised that the variables faob,fnb,fnsp,famx,fhan,SseBF,Sno3BF,Sno2BF,Snh4BF,So2BF are (nx+1)x1 but the results of their products (muaob,munb,munsp,muamx,muhan and muaver) are (nx)x1. How is this possible?

I don't know why MATLAB is not throwing an error in the for loop.

SCRPT:

tstart = 0;  
tend = 10;   
nt = 10;   
xstart = 0;  
xend = 0.012;  
nx = 10;
t = linspace(tstart,tend,nt);   
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,6);
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = 1*10^-6;

y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17]; 
%% execution 
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1);
%writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 =  oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end

FUNCTION:

function dydt=MBBR(t,y)
tstart = 0;  
tend = 10;   
nt = 10;   
xstart = 0;  
xend = 0.012;  
nx = 10;
t = linspace(tstart,tend,nt);   
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; n=0.5; 
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
dydt=@fval;
muaob=u1*(So2BF/(Ko2aob+So2BF))*(Snh4BF./(Knh4aob+Snh4BF));  
MDaobaer =baob*(So2BF./(Ko2aob+So2BF));   MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)))';
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF));                                        MDnbaer= bnb*(So2BF./(Ko2nb+So2BF));      MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF))';
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF));                                     MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF));  MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF))';
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF));      MDamxaer= bamx*(So2BF./(Ko2amx+So2BF));  MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF))';
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)); 
muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF)); 
muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF));                                  MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF))';
muhan=muhanSo2+muhanSno3+muhanSno2;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L*muaver+sigma;
VL=VR-AL*(L + LL); 
z=0.012; zeta=z/L;
U=L*muaver*zeta;
X=5000; %mg/L
if and (ge(t,0),le(t,5))
   Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
 elseif  and (gt(t,5),le(t,9))     
     Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
 elseif and (gt(t,9), le(t,12))
     Snh4o=909;Qo=0.132*VL*10^3/(909*24);
 elseif and (gt(t,12), le(t,15))
     Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
 elseif and (gt(t,15), le(t,21))
     Snh4o=937;Qo=0.185*VL*10^3/(937*24);
 elseif and (gt(t,21),le(t,33))
     Snh4o=937;Qo=0.2286*VL*10^3/(937*24); 
 elseif and (gt(t,33),le(t,35))
     Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
 elseif and (gt(t,35),le(t,37))
     Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
 elseif and (gt(t,37),le(t,49))
     Snh4o=937;Qo=0.348*VL*10^3/(937*24);
 elseif and (gt(t,49),le(t,77))
     Snh4o=937;Qo=0.448*VL*10^3/(937*24);
 elseif and (gt(t,77),le(t,91))
     Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
 elseif and (gt(t,91),le(t,119))
     Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
 elseif and (gt(t,119),le(t,127))
     Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
  elseif and (gt(t,127),le(t,134))
     Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
 elseif and (gt(t,134),le(t,140))
     Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
 elseif and (gt(t,140),le(t,148))
     Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
elseif and (gt(t,148),le(t,155))
     Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
 elseif and (gt(t,155), le(t,161))
     Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
 elseif and (gt(t,161),le(t,169))
     Snh4o=1078;Qo=0.564*VL*10^3/(1078*24);
 else
     Snh4o=928; Qo=0.47*VL*10^3/(928*24);
end
KaL=0;
while So2BL<0.5
    KaL=80;
end
%Si=SseBF,dSno3BF,dSno2BF, dSnh4BF or dSo2BF
%BC-1 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=0; 
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); 
So2BFmin=So2BF(1); 
%SLi=SseBL,dSno3BL,dSno2BL, dSnh4BL,dSo2BL
%BC-2 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=L*DLSse*(SLi-Si)/LL;
SseBFmax=SseBF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL; 
Sno3BFmax=Sno3BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL; 
Sno2BFmax=Sno2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;  
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%% ODEs (1-6&17)
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBFmax);
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BFmax);
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BFmax);
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BFmax);
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BFmax);
dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/fhan(nx))/(KX+(Xs/fhan(nx))))*fhan(nx);
dLdt=L*muaver(nx)*zeta+sigma;
%% PDEs (7-16)
for ix=2:nx
            
    dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(SseBF(ix+1)-SseBF(ix))/L ...
                - n*(1/YNo2han)*muhanSno2(ix).*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix).*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix).*fhan(ix)*X; %Sse consumption
    dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))/L ...
                - n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix).*fhan(ix)*X + (1/1.14)*muamx(ix).*famx(ix)*X + (1/Ynb)*munb(ix).*fnb(ix)*X ...
                + (1/Ynsp)*munsp(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix).*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix).*fnb(ix) ...
                - ((1-fI)/2.86)*MDnspana(ix).*fnsp(ix)*X -  ((1-fI)/2.86)*MDamxana(ix).*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix).*fhan(ix)*X; %Sno3 consumption+generation
    dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))/L ...
                - ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix).*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix).*famx(ix)*X + (1/Yaob)*muaob(ix).*faob(ix)*X ...
                - (1/Ynb)*munb(ix).*fnb(ix)*X - (1/Ynsp)*munb(ix).*fnsp(ix)*X; %Sno2 comsumption+generation
    dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))/L ...
                - (INBM+1/Yaob)*muaob(ix).*faob(ix)*X -(INBM+1/Yamx)*muamx(ix).*famx(ix)*X - INBM*munb(ix).*fnb(ix)*X -INBM*munsp(ix).*fnsp(ix)*X ...
                - INBM*n*muhanSno2(ix).*fhan(ix)*X - INBM*n*muhanSno3(ix).*fhan(ix)*X - INBM*muhanSo2(ix).*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix).*faob(ix)*X ...
                + (INBM-(fI*INXI))*MDnbana(ix).*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix).*famx(ix)*X ...
                + (INBM-(fI*INXI))*MDhanana(ix).*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix).*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix).*fnb(ix)*X ...
                +(INBM-(fI*INXI))*MDnspaer(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix).*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix).*fhan(ix)*X;%Snh4 consumption+generation
    dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L ...
                -((1-Yhaer)/Yhaer)*muhanSo2(ix).*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix).*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix).*fnb(ix)*X ...
                - ((1.14-Ynsp)/Ynsp)*munsp(ix).*fnsp(ix)*X - (1-fI)*MDaobaer(ix).*faob(ix)*X - (1-fI)*MDnbaer(ix).*fnb(ix)*X -(1-fI)*MDnspaer(ix).*fnsp(ix)*X ...
                - (1-fI)*MDamxaer(ix).*famx(ix)*X - (1-fI)*MDhanaer(ix).*fhan(ix);%So2 inflow+consumption
    dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta*UL(ix)).*(faob(ix+1)-faob(ix))/((x(nx+1)-x(nx-1))*L);
    dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta*UL(ix)).*(fnb(ix+1)-fnb(ix))/((x(nx+1)-x(nx-1))*L);
    dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta*UL(ix)).*(fnsp(ix+1)-fnsp(ix))/((x(nx+1)-x(nx-1))*L);
    dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta*UL(ix)).*(famx(ix+1)-famx(ix))/((x(nx+1)-x(nx-1))*L);
    dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta*UL(ix)).*(fhan(ix+1)-fhan(ix))/((x(nx+1)-x(nx-1))*L);
            
end
%%left BC for faob,fnb,fnsp,famx,fhan
    function fval=BC(t,yBC,muaob,munb,munsp,muamx,muhan,muaver)
            fval=zeros(length(yBC),1);
            faobBC=yBC(1,:);
            fnbBC=yBC(2,:);
            fnspBC=yBC(3,:);
            famxBC=yBC(4,:);
            fhanBC=yBC(5,:);
            for i=1
                fval=[(muaob(i)-muaver(i))*faobBC(i)
                      (munb(i)-muaver(i))*fnbBC(i)
                      (munsp(i)-muaver(i))*fnspBC(i)
                      (muamx(i)-muaver(i))*famxBC(i)
                      (muhan(i)-muaver(i))*fhanBC(i)];
            end
            fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
        end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end

[Torsten]

I have relooked at the codes and realised that the variables faob,fnb,fnsp,famx,fhan,SseBF,Sno3BF,Sno2BF,Snh4BF,So2BF are (nx+1)x1 but the results of their products (muaob,munb,munsp,muamx,muhan and muaver) are (nx)x1. How is this possible?

When you read faob,fnb,fnsp,famx,fhan,SseBF,Sno3BF,Sno2BF,Snh4BF,So2BF from the solution vector y, all these arrays become nx x 1 and so are muaob,munb,munsp,muamx,muhan,muaver that you deduce from them. Later on in your code, you most probaby enlarge your solution variable vectors by their boundary value such that they become (nx+1) x 1.

[KIPROTICH KOSGEY]

Thank you Torsten. Your assistance is greatly appreciated, be blessed.

I have taken all your comments into account, though I sometimes come short as I am not a MATLAB expert.

I have moved the codes for computation of muaob,munb,munsp,muamx,muhan and muaver until after the codes for enlargement of faob,fnb,fnsp,famx,fhan,SseBF,Sno3BF,Sno2BF,Snh4BF,So2BF so that all the solution variables become (nx+1)x1. However, this has no effect on the computations since either way MATLAB was self-updating the changes.

As for the endless computations (stretching over 24 hours withouting completing), I have figured out that the while loop caused all these problems:

KaL=0;

while So2BL<0.5

KaL=80;

end

I have then edited these codes as:

KaL=0;

while So2BL<0.5

KaL=80;

break

end

Is this a correct move? Even so, So2BL keeps decreasing yet switching KaL between 0 and 80 is supposed to keep this value within some range.

Also, I am not quite sure on how to incorporate the left BC at x=0 for the pdes d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=0. The integrals are: SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); So2BFmin=So2BF(1);

when I set as SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax]; MATLAB throws an error in the for loop 'index must not exceeed 2' for these variables. So how do i factor in the left BCs?

Sorry to bother this much but I must say that I have learnt a lot from you and Walter throughout this exchange.

Thanks in advance.

[Walter Roberson]

if and (ge(t,0),le(t,5))
   Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
 elseif  and (gt(t,5),le(t,9))     

The mathematics of Runge Kutta ODE solvers requires that the code you provide must have continuous first and second derivatives. When that condition is violated, if you are lucky you get an error message about being unable to integrate at a particular time; if you are not lucky then you are not told that your results are invalid.

Your if/elseif looks to be time dependent. You need to recode as a series of calls to ode45(), one for each block of time that the constants such as Snh4o are consistent for, possibly adjusting the boundary conditions (for example if you were injecting chemicals at particular times.)

[KIPROTICH KOSGEY]

Thank you Torsten and Walter. Despite my efforts I have failed to achieve my aims.

I have fixed Snh4o and Qo so that none of the inputs vary with time but now is now throwing the error: Warning: Failure at t=1.317414e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

I could not also code so that KaL is 80 when So2BL is <0.5 and 0 when above this value. I am giving up. Below is my final program, it might be of help to somebody else.

SCRIPT:

tstart = 0;  
tend = 10*24;   
nt = 30000;   
xstart = 0;  
xend = 0.012;  
nx = 100;
t = linspace(tstart,tend,nt);   
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,6);
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = 1*10^-6;

y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17]; 
%% execution 
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1);

Warning: Failure at t=1.691992e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.440892e-16) at time t.

%writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 =  oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end

Function:

function dydt=MBBR(t,y)
tstart = 0;  
tend = 10;   
nt = 30000;   
xstart = 0;  
xend = 0.012;  
nx = 100;
t = linspace(tstart,tend,nt);   
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; n=0.5; 
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
dydt=@fval;
VL=VR-AL*(L + LL); %m^3
Snh4o=900; %mg/L 
Qo=VL*2; 
KaL=80;
%Si=SseBF,dSno3BF,dSno2BF, dSnh4BF or dSo2BF
%BC-1 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=0; 
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); 
So2BFmin=So2BF(1); 
%SLi=SseBL,dSno3BL,dSno2BL, dSnh4BL,dSo2BL
%BC-2 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=L*DLSse*(SLi-Si)/LL;
SseBFmax=SseBF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL; 
Sno3BFmax=Sno3BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL; 
Sno2BFmax=Sno2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;  
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
muaob=u1*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF));                                   MDaobaer =baob*(So2BF./(Ko2aob+So2BF));   MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)))';
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF));                                        MDnbaer= bnb*(So2BF./(Ko2nb+So2BF));      MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF))';
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF));                                     MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF));  MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF))';
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF));      MDamxaer= bamx*(So2BF./(Ko2amx+So2BF));  MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF))';
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)); 
muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF)); 
muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF));                                  MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF))';
muhan=muhanSo2+muhanSno3+muhanSno2;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L*muaver+sigma;
z=0.012; zeta=z/L;
U=L*muaver*zeta;
X=5000; %mg/L
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%% ODEs (1-6&17)
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBFmax);
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BFmax);
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BFmax);
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BFmax);
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BFmax);
dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/(fhan(nx)*X))/(KX+(Xs/(fhan(nx)*X))))*(fhan(nx)*X);
dLdt=L*muaver(nx)*zeta+sigma;
%% PDEs (7-16)
for ix=2:nx
    
    dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(SseBF(ix+1)-SseBF(ix))/L ...
        - n*(1/YNo2han)*muhanSno2(ix).*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix).*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix).*fhan(ix)*X; %Sse consumption
    dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))/L ...
        - n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix).*fhan(ix)*X + (1/1.14)*muamx(ix).*famx(ix)*X + (1/Ynb)*munb(ix).*fnb(ix)*X ...
        + (1/Ynsp)*munsp(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix).*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix).*fnb(ix) ...
        - ((1-fI)/2.86)*MDnspana(ix).*fnsp(ix)*X -  ((1-fI)/2.86)*MDamxana(ix).*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix).*fhan(ix)*X; %Sno3 consumption+generation
    dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))/L ...
        - ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix).*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix).*famx(ix)*X + (1/Yaob)*muaob(ix).*faob(ix)*X ...
        - (1/Ynb)*munb(ix).*fnb(ix)*X - (1/Ynsp)*munb(ix).*fnsp(ix)*X; %Sno2 comsumption+generation
    dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))/L ...
        - (INBM+1/Yaob)*muaob(ix).*faob(ix)*X -(INBM+1/Yamx)*muamx(ix).*famx(ix)*X - INBM*munb(ix).*fnb(ix)*X -INBM*munsp(ix).*fnsp(ix)*X ...
        - INBM*n*muhanSno2(ix).*fhan(ix)*X - INBM*n*muhanSno3(ix).*fhan(ix)*X - INBM*muhanSo2(ix).*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix).*faob(ix)*X ...
        + (INBM-(fI*INXI))*MDnbana(ix).*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix).*famx(ix)*X ...
        + (INBM-(fI*INXI))*MDhanana(ix).*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix).*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix).*fnb(ix)*X ...
        +(INBM-(fI*INXI))*MDnspaer(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix).*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix).*fhan(ix)*X;%Snh4 consumption+generation
    dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L ...
        -((1-Yhaer)/Yhaer)*muhanSo2(ix).*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix).*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix).*fnb(ix)*X ...
        - ((1.14-Ynsp)/Ynsp)*munsp(ix).*fnsp(ix)*X - (1-fI)*MDaobaer(ix).*faob(ix)*X - (1-fI)*MDnbaer(ix).*fnb(ix)*X -(1-fI)*MDnspaer(ix).*fnsp(ix)*X ...
        - (1-fI)*MDamxaer(ix).*famx(ix)*X - (1-fI)*MDhanaer(ix).*fhan(ix);%So2 inflow+consumption
    dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta*UL(ix)).*(faob(ix+1)-faob(ix))/((x(nx+1)-x(nx-1))*L);
    dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta*UL(ix)).*(fnb(ix+1)-fnb(ix))/((x(nx+1)-x(nx-1))*L);
    dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta*UL(ix)).*(fnsp(ix+1)-fnsp(ix))/((x(nx+1)-x(nx-1))*L);
    dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta*UL(ix)).*(famx(ix+1)-famx(ix))/((x(nx+1)-x(nx-1))*L);
    dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta*UL(ix)).*(fhan(ix+1)-fhan(ix))/((x(nx+1)-x(nx-1))*L);
    
end
%%left BC for faob,fnb,fnsp,famx,fhan
    function fval=BC(t,yBC,muaob,munb,munsp,muamx,muhan,muaver)
        fval=zeros(length(yBC),1);
        faobBC=yBC(1,:);
        fnbBC=yBC(2,:);
        fnspBC=yBC(3,:);
        famxBC=yBC(4,:);
        fhanBC=yBC(5,:);
        for i=1
            fval=[(muaob(i)-muaver(i))*faobBC(i)
                (munb(i)-muaver(i))*fnbBC(i)
                (munsp(i)-muaver(i))*fnspBC(i)
                (muamx(i)-muaver(i))*famxBC(i)
                (muhan(i)-muaver(i))*fhanBC(i)];
        end
        fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
    end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end

Thanks again and be blessed!

[Walter Roberson]

If you switch to ode23s then you can get up to about 1.32 seconds.

Once you do that, get busy plotting. I suggest plotting 10 at a time, such as

NN = size(ySol,2);
for K = 1 : 10 : NN
    idx = K:min(K+9, NN);
    figure
    plot(tSol, ySol(:,idx));
    legend("y_{" + idx + "}");
end

You will see that starting roughly ySol(:,510) that every 100 goes quite steep -- e.g., 519, 619, 719-ish are all steep.

But before that, starting from the second group of 10, ySol(:,11:20) that every plot goes steep around 1.3 for at least some of the outputs.

I am not going to analyze the code to figure out "why" . When it comes to ODEs, sometimes the reason is just "because it does" -- that either the equations are unstable (common!) or else that the equations need extended precision calculations to remain stable (certainly happens, but a lot of the time that loss of precision is suspected as the cause, it turns out that the equations are unstable and high precision calculations just postpones the time thte problem shows up at.)

[KIPROTICH KOSGEY]

I am trying to simulate a biofilm process in which substrates diffuse through some thick layers where microorganisms are residing. This means that none of the variables should be negative. In addition, the key substrate So2 need to range between 0.5 and 1.5 mg/L. To control this value, I have was initially trying to switch KaL between 0 and 80: when So2BL<0.5, KaL=80, when So2BL>1.5, KaL=0, hence the while loop in the below codes. However, MATLAB wzs ignoring this for some reason:

while So2BL<0.5

KaL=80;

break

end

I also included 'NonNegative',1 in options so that none of the solution variables could be negative but then this also didn't help.

I therefore feel that if I could succeeed in setting the two parameters then I could stabilise the calculations.

[Walter Roberson]

Use Event functions to at least detect the case where the values go out of physical range. But you will need to decide whether you should just terminate the system there (that is, it might be telling you that the equations or initial conditions must be wrong) — or if instead you want to adjust the system (for example corresponding to adding chemicals or changing temperature or the like.)

[KIPROTICH KOSGEY]

Dear Walter,

Your comments and suggestions are highly appreciated.

I have gone through the codes and didn't find any error relating to the equations. However, I see that faob goes to negative leading to unstable calculations relating to So2BF as a result of -*(-) leading to + when in fact the it should be -.

As suggested, I have introduced an event function in the script but MATLAB doesn't execute it, possibly because of my horrible coding style. I would like to stop the execution when any of the solution variables go to zero.

Sript:

t=[0 180*24];
t1=[0 5*24 9*24 12*24 15*24 21*24 33*24 35*24 37*24 49*24 77*24 91*24 119*24 127*24 134*24 140*24 148*24 155*24 161*24 169*24 180*24];
xstart = 0;  
xend = 0.012;
nx=100;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
nx=100;
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5);
y06 = b(6,:).'; y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).';
y12 = b(12,:).'; y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = 1*10^-6;

y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16]; 
%% execution 
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1,options);
%writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 =  oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
function [value,isterminal,direction]=eventsfcn(t,y)
value=[y(1);y(2);y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16)];
isterminal = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];        % stop at 0
direction  = [-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1];    % decreasing
end

Function

function dydt=MBBR(t,y)
nx=100;
xstart = 0;  
xend = 0.012;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];

%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han 
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;  
YNo3han=0.5*0.79; Yhaer=0.37; 
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; n=0.5; 
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5;  DLSno3=1.8*10^-4;  DLSno2=1.8*10^-4;  DLSnh4=1.8*10^-4;  DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5;   DSno3=1.4*10^-5;   DSno2=1.4*10^-4;   DSnh4=1.5*10^-4;   DSo2=2.2*10^-4; 
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
SseBF=y(6:nx+5,:);
Sno3BF=y(nx+6:2*nx+5,:);
Sno2BF=y(2*nx+6:3*nx+5,:);
Snh4BF=y(3*nx+6:4*nx+5,:);
So2BF=y(4*nx+6:5*nx+5,:);
faob=y(5*nx+6:6*nx+5,:);
fnb=y(6*nx+6:7*nx+5,:);
fnsp=y(7*nx+6:8*nx+5,:);
famx=y(8*nx+6:9*nx+5,:);
fhan=y(9*nx+6:10*nx+5,:);
L=y(10*nx+6,:);
VL=VR-AL*(L + LL); %m^3
Snh4o=900; %mg/L 
Qo=VL*2; 
KaL=80;
%Si=SseBF,dSno3BF,dSno2BF, dSnh4BF or dSo2BF
%BC-1 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=0; 
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); 
So2BFmin=So2BF(1); 
%SLi=SseBL,dSno3BL,dSno2BL, dSnh4BL,dSo2BL
%BC-2 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=L*DLSse*(SLi-Si)/LL;
SseBFmax=SseBF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL; 
Sno3BFmax=Sno3BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL; 
Sno2BFmax=Sno2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL; 
Snh4BFmax=Snh4BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;  
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
muaob=u1*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF));                                   MDaobaer =baob*(So2BF./(Ko2aob+So2BF));   MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)))';
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF));                                        MDnbaer= bnb*(So2BF./(Ko2nb+So2BF));      MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF))';
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF));                                     MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF));  MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF))';
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF));      MDamxaer= bamx*(So2BF./(Ko2amx+So2BF));  MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF))';
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)); 
muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF)); 
muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF));                                  MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF))';
muhan=muhanSo2+muhanSno3+muhanSno2;
kinert=0.1; %day^-1
finert=1-(faob+fnb+fnsp+famx+fhan);
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan+finert.*kinert*5;
UL=L*muaver+sigma;
z=0.012; zeta=z/L;
U=L*muaver*zeta;
X=5000; %mg/L 

if and (ge(t,0),le(t,5))
    Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
elseif  and (gt(t,5),le(t,9))     
    Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
elseif and (gt(t,9), le(t,12))
    Snh4o=909;Qo=0.132*VL*10^3/(909*24);
elseif and (gt(t,12), le(t,15))
    Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
elseif and (gt(t,15), le(t,21))
    Snh4o=937;Qo=0.185*VL*10^3/(937*24);
elseif and (gt(t,21),le(t,33))
    Snh4o=937;Qo=0.2286*VL*10^3/(937*24); 
elseif and (gt(t,33),le(t,35))
    Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
elseif and (gt(t,35),le(t,37))
    Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
elseif and (gt(t,37),le(t,49))
    Snh4o=937;Qo=0.348*VL*10^3/(937*24);
elseif and (gt(t,49),le(t,77))
    Snh4o=937;Qo=0.448*VL*10^3/(937*24);
elseif and (gt(t,77),le(t,91))
    Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
elseif and (gt(t,91),le(t,119))
    Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
elseif and (gt(t,119),le(t,127))
    Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
elseif and (gt(t,127),le(t,134))
    Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
elseif and (gt(t,134),le(t,140))
    Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
elseif and (gt(t,140),le(t,148))
    Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
elseif and (gt(t,148),le(t,155))
    Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
elseif and (gt(t,155), le(t,161))
    Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
elseif and (gt(t,161),le(t,169))
    Snh4o=1078;Qo=0.564*VL*10^3/(1078*24);
else
    Snh4o=928; Qo=0.47*VL*10^3/(928*24);
end
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%% ODEs (1-6&17)
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BF(nx));
dLdt=L*muaver(nx)*zeta+sigma;
%% PDEs (7-16)
for ix=2:nx
    
    dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(SseBF(ix+1)-SseBF(ix))/L ...
        - n*(1/YNo2han)*muhanSno2(ix).*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix).*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix).*fhan(ix)*X; %Sse consumption
    dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))/L ...
        - n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix).*fhan(ix)*X + (1/1.14)*muamx(ix).*famx(ix)*X + (1/Ynb)*munb(ix).*fnb(ix)*X ...
        + (1/Ynsp)*munsp(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix).*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix).*fnb(ix) ...
        - ((1-fI)/2.86)*MDnspana(ix).*fnsp(ix)*X -  ((1-fI)/2.86)*MDamxana(ix).*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix).*fhan(ix)*X; %Sno3 consumption+generation
    dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))/L ...
        - ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix).*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix).*famx(ix)*X + (1/Yaob)*muaob(ix).*faob(ix)*X ...
        - (1/Ynb)*munb(ix).*fnb(ix)*X - (1/Ynsp)*munb(ix).*fnsp(ix)*X; %Sno2 comsumption+generation
    dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))/L ...
        - (INBM+1/Yaob)*muaob(ix).*faob(ix)*X -(INBM+1/Yamx)*muamx(ix).*famx(ix)*X - INBM*munb(ix).*fnb(ix)*X -INBM*munsp(ix).*fnsp(ix)*X ...
        - INBM*n*muhanSno2(ix).*fhan(ix)*X - INBM*n*muhanSno3(ix).*fhan(ix)*X - INBM*muhanSo2(ix).*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix).*faob(ix)*X ...
        + (INBM-(fI*INXI))*MDnbana(ix).*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix).*famx(ix)*X ...
        + (INBM-(fI*INXI))*MDhanana(ix).*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix).*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix).*fnb(ix)*X ...
        +(INBM-(fI*INXI))*MDnspaer(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix).*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix).*fhan(ix)*X;%Snh4 consumption+generation
    dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L ...
        -((1-Yhaer)/Yhaer)*muhanSo2(ix).*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix).*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix).*fnb(ix)*X ...
        - ((1.14-Ynsp)/Ynsp)*munsp(ix).*fnsp(ix)*X - (1-fI)*MDaobaer(ix).*faob(ix)*X - (1-fI)*MDnbaer(ix).*fnb(ix)*X -(1-fI)*MDnspaer(ix).*fnsp(ix)*X ...
        - (1-fI)*MDamxaer(ix).*famx(ix)*X - (1-fI)*MDhanaer(ix).*fhan(ix);%So2 inflow+consumption
    dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta*UL(ix)).*(faob(ix+1)-faob(ix))/((x(nx+1)-x(nx-1))*L);
    dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta*UL(ix)).*(fnb(ix+1)-fnb(ix))/((x(nx+1)-x(nx-1))*L);
    dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta*UL(ix)).*(fnsp(ix+1)-fnsp(ix))/((x(nx+1)-x(nx-1))*L);
    dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta*UL(ix)).*(famx(ix+1)-famx(ix))/((x(nx+1)-x(nx-1))*L);
    dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta*UL(ix)).*(fhan(ix+1)-fhan(ix))/((x(nx+1)-x(nx-1))*L);
    
end
%%left BC for faob,fnb,fnsp,famx,fhan
dydt=@fval;
    function fval=BC(t,yBC,muaob,munb,munsp,muamx,muhan,muaver)
        fval=zeros(length(yBC),1);
        faobBC=yBC(1,:);
        fnbBC=yBC(2,:);
        fnspBC=yBC(3,:);
        famxBC=yBC(4,:);
        fhanBC=yBC(5,:);
        for i=1
            fval=[(muaob(i)-muaver(i))*faobBC(i)
                (munb(i)-muaver(i))*fnbBC(i)
                (munsp(i)-muaver(i))*fnspBC(i)
                (muamx(i)-muaver(i))*famxBC(i)
                (muhan(i)-muaver(i))*fhanBC(i)];
        end
        fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
    end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end

[Torsten]

y06-y15 all have nx elements, not only one.

Thus

value=[y(1);y(2);y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16)];
isterminal = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];        % stop at 0
direction  = [-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1];    % decreasing

is wrong.

But how should the event function help in your case ? The solver will stop integrating almost instantly. If you don't have a strategy on how to continue after the event occurs, you are left with a solution vector after 1 sec or so.

[KIPROTICH KOSGEY]

Dear Torsten and Walter.

Thank you so much for giving me hints on possible solutions to the coding mistakes I have been making. I have learnt a lot in the course of this exchange.

I have since edited the script as shown below and coded for the program when So2BL is above or below the given range (since this variable affects all the other variables directly leading to unstable calculations and failures). I have also divided the functions into two to gather for different periods (aerobic (increasing So2BL because KaL=80) and anaerobic (decreasing So2BL because KaL=0)). Although it apperas to be working, I keep getting the error: out of memory. Could there be a way to go around this?

Main functions attached.

Script:

tend=180*24;   
t=[0 180*24];
t1=[0 5*24 9*24 12*24 15*24 21*24 33*24 35*24 37*24 49*24 77*24 91*24 119*24 127*24 134*24 140*24 148*24 155*24 161*24 169*24 180*24];
xstart = 0;  
xend = 0.012;
nx=100;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
nx=100;
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5); y16 = 1*10^-6;
y06 = b(6,:).'; y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).';
y12 = b(12,:).'; y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; 
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16]; 
%% execution 
options1 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_aerobic);
options2 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_anoxic);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:2
    y2=[ySol(end,1);ySol(end,2);ySol(end,3);ySol(end,4);ySol(end,5);ySol(end,6:105)';ySol(end,106:205)';ySol(end,206:305)';ySol(end,306:405)';...
        ySol(end,406:505)';ySol(end,506:605)';ySol(end,606:705)';ySol(end,706:805)';ySol(end,806:905)';ySol(end,906:1005)';ySol(end,1006)'];
    [tSol,ySol]=ode45(@MBBR_anoxic,[tplot(end) tend],y2,options2);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
z1=ceil(1.5-z);
value=z1;
isterminal = 1;    % stop at 0
direction  = -1;    % increasing
end
%% events anoxic phase
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
z=y(5);
z1=ceil(0.5-z);
value=z1;
isterminal = 1;    % stop at 0
direction  = 1;    % decreasing
end

[Torsten]

In which phase of execution of your program do you get the error message ? While using ode15s or ode45 ? If ode45 is called with a two-element vector for tspan as you do, it will collect so many output solutions that you will get a memory overflow.

Why do you use ode45 at all ? It is not suited for your kind of problem.

[KIPROTICH KOSGEY]

I was getting an error in the second phase of execution and I agree with you that ode45 was causing this error! I don't even know why I used ode45 in the first plast. I should have also used yplot instead of ySol in the initialisation of the second phase.

However, after correcting this, I have found that the program is not calling the second function MBBR_anoxic which is supposed to lead to a decrease in So2BL, but rather continues with the first program MBBR_aerobic leading to an exponential increasing in So2BL until the program fails.

[Torsten]

However, after correcting this, I have found that the program is not calling the second function MBBR_anoxic which is supposed to lead to a decrease in So2BL, but rather continues with the first program MBBR_aerobic leading to an exponential increasing in So2BL until the program fails.

Seems you missed to switch a + to a - somewhere in your code for the second phase.

[KIPROTICH KOSGEY]

Thank you Torsten, I found the problem leading to the below issue. However, after editting it and modifying the script, MATLAB stops just below y(5) is 1.5 yet I want it to stop when y(5)>1.5.

t=[0 180*24];
xstart = 0;  
xend = 0.012;
nx=100;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
nx=100;
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5); y16 = 1*10^-6;
y06 = b(6,:).'; y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).';
y12 = b(12,:).'; y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; 
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16]; 
%% execution 
options1 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_aerobic);
options2 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_anoxic);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:2
    y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
        yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
    [tSol,ySol,tev,yev]=ode15s(@MBBR_anoxic,[tplot(end) 5*180],y2,options2);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
z1=+(z<=1.5);
value=z1;
isterminal = 1;    % stop at 0
direction  = [];    % decreasing
end
%% events anoxic phase
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
e=y(5);
e1=+(e>0.5);
value=e1;
isterminal = 1;    % stop at 0
direction  = [];    % decreasing
end

[KIPROTICH KOSGEY]

Never mind, I found the cause of the failure. Thanks again!

[fahad isah adam]

Kudos to you, Sir.

[Torsten]

The event function you define must be differentiable with respect to the solution variables. z1 and e1 are not differentiable, but jump functions.

[KIPROTICH KOSGEY]

I am learning as I have always done from the time I sought assistance regarding my codes.

with z1 and e1, the program is running and the variable y(5) is changing as expected. However, on the third round, the solution variables are curiously constant when they should have been changing. Could this be related to the issue you have raised?

Basically, I would like to keep y(5) between 0.5 and 1.5. An event as per the examples in MATLAB documentation is recorded when a function crosses zero. So in my case, how can I meet this requirement whereas my variables range above this?

This is my script:

close all; clear; clc;
t1=[0 180*24];
xstart = 0;  
xend = 0.012;
nx=100;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
nx=100;
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5); y16 = 1*10^-6;
y06 = b(6,:).'; y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).';
y12 = b(12,:).'; y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; 
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16]; 
%% execution 
options1 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_aerobic);
options2 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_anoxic);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t1,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:20
    y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
        yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
    [tSol,ySol]=ode15s(@MBBR_anoxic,[tplot(end) 24*180],y2,options2);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
    y3=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
        yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
    [tSol,ySol]=ode15s(@MBBR_aerobic,[tplot(end) 24*180],y3,options1);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
z1=+(z<=1.5);
value=z1;
isterminal = 1;    % stop at 0
direction  = -1;    % decreasing
end
%% events anoxic phase
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
e=y(5);
e1=+(e>0.5);
value=e1;
isterminal = 1;    % stop at 0
direction  = -1;    % decreasing
end 

[Walter Roberson]

e1=+(e>0.5);

That is not differentiatable. Use

value = [0.5 - e, e - 1.5];

MATLAB will attempt to keep both of those non-positive (so <= 0) .

0.5 - e is less than or equal to 0 when e is 0.5 or greater.

e - 1.5 is less than or equal to 0 when e is no more than 1.5.

Therefore the [] I show will confine e to be in the range [0.5, 1.5]

[KIPROTICH KOSGEY]

Thank you Walter and Torsten.

I have edited my script and the program is now running.

However, Since I am trying to call a different function each time an event occurs, is it the correct coding to do the below. I am doing this to mimic process controller that regulates aeration in a biological system that switches aeration on and off depending on the concentration of oxygen y(5) in the bulk solution. The difference between the two functions is that KaL is either 0 or 10 in MBBR-aerobic and MBBR-anoxic:

close all; clear; clc;
t1=[0 180*24];
xstart = 0;  
xend = 0.012;
nx=100;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
nx=100;
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5); y16 = 1*10^-6;
y06 = b(6,:).'; y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).';
y12 = b(12,:).'; y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; 
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16]; 
%% execution 
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t1,y1,options);
tplot=tSol;
yplot=ySol;
for i=1:200
    y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
        yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
    [tSol,ySol]=ode15s(@MBBR_anoxic,[tplot(end) 24*180],y2,options);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
    y3=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
        yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
    [tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),[tplot(end) 24*180],y3,options);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn(t,y)
z=y(5);
value=[0.5 - z; z - 1.5];
isterminal = [1;1] ;   % stop at 0
direction  = [1;1];    % increasing
end

[KIPROTICH KOSGEY]

Dear Torsten and Walter,

Thank you so much for the several comments and insights into coding with Matlab.

I have improved my code further by adding an integral for UL and U but its giving me the error:

Error using integral Limits of integration must be double or single scalars.

This are the lines of codes:

UL=integral(@(L) muaver,0,L,'ArrayValued',true)+sigma;

U=integral(@(x) muaver,0,x,'ArrayValued',true);

[KIPROTICH KOSGEY]

Oh sorry, I have realised the mistake I made, I used x instead of xend=0.012 in the limits of integration.

sorry again to have bothered you.

Regards

[KIPROTICH KOSGEY]

Hi,

This is my current script. I have set 'NoNegative', 1, but still the y values are going to negative, what could be the problem?

close all; clear; clc;
dbstop if error
t1=[0 180*24];
xstart = 0; 
xend = 0.012;
nx=10;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5); y16 = 1*10^-6;
y06 = b(6,:).'; y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).';
y12 = b(12,:).'; y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; 
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16]; 
%% execution 
options1 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_aerobic,'NonNegative',1);
options2 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_anoxic,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t1,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:10
    y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:15)';yplot(end,16:25)';yplot(end,26:35)';yplot(end,36:45)';...
        yplot(end,46:55)';yplot(end,56:65)';yplot(end,66:75)';yplot(end,76:85)';yplot(end,86:95)';yplot(end,96:105)';yplot(end,106)];
    [tSol,ySol]=ode15s(@MBBR_anoxic,[tplot(end) 24*180],y2,options2);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
    y3=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:15)';yplot(end,16:25)';yplot(end,26:35)';yplot(end,36:45)';...
        yplot(end,46:55)';yplot(end,56:65)';yplot(end,66:75)';yplot(end,76:85)';yplot(end,86:95)';yplot(end,96:105)';yplot(end,106)];
    [tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),[tplot(end) 24*180],y3,options1);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [10;10;10;10;1;10;10;10;10;10;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
value=z - 1.5;
isterminal =1;   % stop at 0
direction  = 1;    % increasing
end
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
z=y(5);
value=0.5 - z;
isterminal = 1 ;   % stop at 0
direction  = 1;    % increasing
end

[Torsten]

If your differential equations are such that the solution becomes negative, you cannot prevent this with the NonNegative option. E.g. the solution of the differential equation y' = -1 with y(0) = 1 is y(t) = 1 - t, and nothing can stop y(t) to become negative for t > 1. It will only be helpful to prevent solution variables to become negative because of numerical issues.

This having set, 'NonNegative',1 only sets the NonNegative option for y(1). To set it for certain or all components of the vector of solution variables, you have to list them all instead of only the 1.

[KIPROTICH KOSGEY]

oh thank you Torsten.

my system consists of odes (y(1:5&16) and pdes (y(6:15)). Would it be correct to set 'NonNegative', [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;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] ?

Or how best to take care of the solution variables for pdes in the mesh?

[Torsten]

odeset('NonNegative',1:6+10*nx)

if you want to keep all your solution components positive.

But as I said: I think your equations are such that the solution becomes negative from your setting - NonNegative cannot help you in this case.

[Walter Roberson]

The 'Nonnegative' option is, as @Torsten indicates, only intended to compensate for numeric roundoff. It does not properly reset boundary conditions to ensure that the integration is valid. It does not have the effect of "bouncing off of 0" (like a ball might bounce on a floor.)

[Walter Roberson]

Is there a newer version of your other .m files since https://www.mathworks.com/matlabcentral/answers/2007287-index-exceeds-the-number-of-array-elements-index-must-not-exceed-1#comment_2866821 ?

Index in position 1 exceeds array bounds. Index must not exceed 106.

Error in MBBR_aerobic (line 49)

Sno3BF=y(nx+6:2*nx+5,:);

y is length 106, nx is 100 so nx+6:2*nx+5 would be 106:205 which is mostly past the end of y

[KIPROTICH KOSGEY]

Dear Walter and Torsten,

Thank you so much for the useful comments.

Including odeset('NonNegative',1:6+10*nx) has prevented the solutions from going to negative but the computation keeps failing. Removing this conditions is not helping either as the computation still fail before it reaches tend.

Either way I am not winning.

This is the script and the main files are atatched:

close all; clear; clc;
dbstop if error
t1=[0 180*24];
xstart = 0; 
xend = 0.012;
nx=10;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5); y16 = 1*10^-6;
y06 = b(6,:).'; y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).';
y12 = b(12,:).'; y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; 
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16]; 
%% execution 
options1 = odeset('RelTol',1e-14,'AbsTol',1e-16,'events',@eventsfcn_aerobic,'NonNegative',1:6+10*nx);
options2 = odeset('RelTol',1e-14,'AbsTol',1e-16,'events',@eventsfcn_anoxic,'NonNegative',1:6+10*nx);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t1,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:10
    y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:15)';yplot(end,16:25)';yplot(end,26:35)';yplot(end,36:45)';...
        yplot(end,46:55)';yplot(end,56:65)';yplot(end,66:75)';yplot(end,76:85)';yplot(end,86:95)';yplot(end,96:105)';yplot(end,106)];
    [tSol,ySol]=ode15s(@MBBR_anoxic,[tplot(end) 24*180],y2,options2);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
    y3=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:15)';yplot(end,16:25)';yplot(end,26:35)';yplot(end,36:45)';...
        yplot(end,46:55)';yplot(end,56:65)';yplot(end,66:75)';yplot(end,76:85)';yplot(end,86:95)';yplot(end,96:105)';yplot(end,106)];
    [tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),[tplot(end) 24*180],y3,options1);
    tplot=[tplot;tSol];
    yplot=[yplot;ySol];
end
%% initial conditions
function u0 =  oscic(x)
u0 = [10;10;10;10;0.5;10;10;10;10;10;0.5;0.4;0.05;0.05;0.3;0.2;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
value=z - 1.5;
isterminal =1;   % stop at 0
direction  = 1;    % increasing
end
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
z=y(5);
value=0.5 - z;
isterminal = 1 ;   % stop at 0
direction  = 1;    % increasing
end