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

Question

KIPROTICH KOSGEY 2023-8-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2007287-index-exceeds-the-number-of-array-elements-index-must-not-exceed-1

编辑： Walter Roberson 2023-9-5

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

52 个评论
显示 50更早的评论隐藏 50更早的评论

KIPROTICH KOSGEY 2023-8-12

编辑：Walter Roberson 2023-8-16

在 MATLAB Online 中打开

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

KIPROTICH KOSGEY 2023-8-15

编辑：Walter Roberson 2023-8-16

在 MATLAB Online 中打开

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 2023-8-16

编辑：Walter Roberson 2023-8-16

在 MATLAB Online 中打开

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 2023-8-16

编辑：Walter Roberson 2023-8-16

在 MATLAB Online 中打开

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 2023-8-25

编辑：Torsten 2023-8-25

在 MATLAB Online 中打开

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

KIPROTICH KOSGEY 2023-8-25

编辑：Torsten 2023-8-25

在 MATLAB Online 中打开

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

Torsten 2023-8-25

编辑：Torsten 2023-8-25

在 MATLAB Online 中打开

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 2023-8-28

编辑：Torsten 2023-8-28

在 MATLAB Online 中打开

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

KIPROTICH KOSGEY 2023-8-28

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.

KIPROTICH KOSGEY 2023-8-28

编辑：Walter Roberson 2023-8-28

在 MATLAB Online 中打开

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 2023-8-28

在 MATLAB Online 中打开

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 2023-8-30

编辑：Walter Roberson 2023-8-30

在 MATLAB Online 中打开

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

KIPROTICH KOSGEY 2023-8-30

编辑：Walter Roberson 2023-8-30

在 MATLAB Online 中打开

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

KIPROTICH KOSGEY 2023-9-1

编辑：Walter Roberson 2023-9-1

在 MATLAB Online 中打开

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 2023-9-1

编辑：Walter Roberson 2023-9-1

在 MATLAB Online 中打开

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 

KIPROTICH KOSGEY 2023-9-2

编辑：Walter Roberson 2023-9-2

在 MATLAB Online 中打开

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 2023-9-4

编辑：Walter Roberson 2023-9-4

在 MATLAB Online 中打开

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

Walter Roberson 2023-9-4

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 2023-9-5

编辑：Walter Roberson 2023-9-5

在 MATLAB Online 中打开

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