t = linspace(tstart,tend,nt);
x = linspace(xstart,xend,nx);
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
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];
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;
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;
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;
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;
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;
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;
Sseo=300; Sno3o=0.2; Xso=355; Sno2o=0.2; So2o=0;
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;
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;
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);
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)
dydt = zeros(length(y),1);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
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;
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);
SseBFmax=SseBL; Sno3BFmax=Sno3BL; Sno2BFmax=Sno2BL; Snh4BFmax=Snh4BL;
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];
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);
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;
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
u0 = [1;1;1;1;1;1;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];