How do I fix 'Error using odearguments @.... returns a vector of length 0, but the length of initial conditions vector is 24...' error?

Hi there,

I keep getting this error and I can't find where the problem is. Any hint would be appreciated.

Actual error: error using odearguments @(T,Y)MBBRPDE(T,Y) returns a vector of length 0, but the length of initial conditions vector is 24. The vector returned by @(T,Y)MBBRPDE(T,Y) and the initial conditions vector must have the same number of elements.

%SCRIPT:
%initial conditions
y0=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0.4;0.2;0.1;0.3;0.1;1;1;1]';
tSpan=[0 180*24];
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'InitialStep',1e-10, 'MStateDependence','strong');
[tSol, ySol]=ode15s(@(t,y) MBBRPDE(t,y), tSpan, y0, options);

Error using odearguments
@(T,Y)MBBRPDE(T,Y) returns a vector of length 0, but the length of initial conditions vector is 24. The vector returned by @(T,Y)MBBRPDE(T,Y) and the initial conditions vector must have the same number of elements.

Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

%Function:
function fval=MBBRPDE(t,y)
%Define the three variables
SseBL=y(1);
Sno3BL=y(2);
Sno2BL=y(3);
Snh4BL=y(4);
So2BL=y(5);
Xs=y(6);
muaob=y(7);
munb=y(8);
munsp=y(9);
muamx=y(10);
muhan=y(11);
SseBF=y(12);
Sno3BF=y(13);
Sno2BF=y(14);
Snh4BF=y(15);
So2BF=y(16);
faob=y(17);
fnb=y(18);
fnsp=y(19);
famx=y(20);
fhan=y(21);
muaver=y(22);
UL=y(23);
VL=y(24);
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;
z1=10^-6;
N=length(UL);
hz=N/z1;

SseBL=zeros(N,1); Sno3BL=zeros(N,1); Sno2BL=zeros(N,1); Snh4BL=zeros(N,1); So2BL=zeros(N,1); Xs=zeros(N,1);
muaob=zeros(N,1); munb=zeros(N,1); munsp=zeros(N,1);  muamx=zeros(N,1); muhan=zeros(N,1); SseBF=zeros(N,1);
Sno3BF=zeros(N,1); Sno2BF=zeros(N,1); Snh4BF=zeros(N,1); So2BF=zeros(N,1); faob=zeros(N,1); fnb=zeros(N,1);
fnsp=zeros(N,1); famx=zeros(N,1); fhan=zeros(N,1); muaver=zeros(N,1); UL=zeros(N,1); VL=zeros(N,1);
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;
if t>=0 && t<5
    Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
elseif  t>=5 && t<9       
    Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
elseif t>=9 && t<12
    Snh4o=909;Qo=0.132*VL*10^3/(909*24);
elseif t>=12 && t<=15
    Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
elseif t>15 && t<=21
    Snh4o=937;Qo=0.185*VL*10^3/(937*24);
elseif t>21 && t<=33
    Snh4o=937;Qo=0.2286*VL*10^3/(937*24); 
elseif t>33 && t<=35
    Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
elseif t>35 && t<=37
    Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
elseif t>37 && t<=49
    Snh4o=937;Qo=0.348*VL*10^3/(937*24);
elseif t>49 && t<=77
    Snh4o=937;Qo=0.448*VL*10^3/(937*24);
elseif t>77 && t<91
    Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
elseif t>=91 && t<=119
    Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
elseif t>119 && t<=127
    Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
elseif t>127 && t<=134
    Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
elseif t>134 && t<=140
    Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
elseif t>140 && t<=148
    Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
elseif t>148 && t<=155
    Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
elseif t>155 && t<=161
    Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
elseif t>161 && 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=80;
AL=120400; LL=1;
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; 

%% PDEs
for i=2:hz:N-1
    %BC-1
    SseBF(1)=0; Sno3BF(1)=0; Sno2BF(1)=0; Snh4BF(1)=0; So2BF(1)=0; faob(1)=0.4;fnb(1)=0.1;fnsp(1)=0.1;famx(1)=0.3;fhan(1)=0.1;
    %BC-2
    SseBF(end)=L*DLSse.*(SseBL(end)-SseBF(end))./LL; Sno3BF(end)=L*DLSno3.*(Sno3BL(end)-Sno3BF(end))/LL; Sno2BF(end)=L*DLSno2.*(Sno2BL(end)-Sno2BF(end))/LL; Snh4BF(end)=L*DLSnh4.*(Snh4BL(end)-Snh4BF(end))/LL;...
        So2BF(end)=L*DLSo2.*(So2BL(end)-So2BF(end))./LL; SseBF(end)=SseBL(end); Sno3BF(end)=Sno3BL(end); Sno2BF(end)=Sno2BL(end); Snh4BF(end)=Snh4BL(end); So2BF(end)=So2BL(end);
    SseBL(i)=Qo.*(Sseo-SseBL(i))./VL(i)-AL.*((DSse/LL)-UL(i)).*(SseBL(i)-SseBF(i));
    Sno3BL(i)=Qo.*(Sno3o-Sno3BL(i))/VL(i)-AL.*((DSno3/LL)-UL(i)).*(Sno3BL(i)-Sno3BF(i));
    Sno2BL(i)=Qo.*(Sno2o-Sno2BL(i))/VL(i)-AL.*((DSno2/LL)-UL(i)).*(Sno2BL(i)-Sno2BF(i));
    Snh4BL(i)=Qo.*(Snh4o-Snh4BL(i))/VL(i)-AL.*((DSnh4/LL)-UL(i)).*(Snh4BL(i)-Snh4BF(i));
    So2BL(i)=Qo.*(So2o-So2BL(i))/VL(i) + KaL.*(8-So2BL(i)) - AL.*((DSo2)-UL(i)).*(So2BL(i)-So2BF(i));
    Xs(i)=Qo.*(Xso-Xs(i))/VL(i) - KH.*((Xs(i)./fhan(i))./(KX+(Xs(i)/fhan(i)))).*fhan(i);
    muaob(i)=u1.*(So2BF(i)./(Ko2aob+So2BF(i))).*(Snh4BF(i)./(Knh4aob+Snh4BF(i))).*faob(i) - baob.*(So2BF(i)./(Ko2aob+So2BF(i))).*faob(i) - baob.*naob.*((Ko2aob./(Ko2aob+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3aob+Sno2BF(i)+Sno3BF(i))).*faob(i);
    munb(i)=u2.*(Sno2BF(i)./(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i) - bnb.*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i) - bnb*nnb.*((Ko2nb./(Ko2nb+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nb+Sno2BF(i)+Sno3BF(i))).*fnb(i);
    munsp(i)=u3.*(Sno2BF(i)./(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i) - bnsp.*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i) - bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nsp+Sno2BF(i)+Sno3BF(i))).*fnsp(i);
    muamx(i)=u5.*((Snh4BF(i)./(Knh4amx+Snh4BF(i))).*(Sno2BF(i)./(Kno2amx+Sno2BF(i))).*(Ko2amx/(Ko2amx+So2BF(i)))).*famx(i) - bamx.*(So2BF(i)./(Ko2amx+So2BF(i))).*famx(i) - bamx.*namx.*((Ko2amx./(Ko2amx+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3amx+Sno2BF(i)+Sno3BF(i))).*famx(i);
    muhan(i)=(u6.*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno2BF(i)./(Kno2han+Sno2BF(i))).*(Ko2han/(Ko2han+So2BF(i))) + u6.*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno3BF(i)./(Kno3han+Sno3BF(i))).*(Ko2han./(Ko2han+So2BF(i)))+u6.*((SseBF(i)./(KSsehan+SseBF(i))).*(So2BF(i)./(Ko2han+So2BF(i)))).*fhan(i))- bhan.*(So2BF(i)./(Ko2han+So2BF(i))).*fhan(i) - bhan.*nhan.*((Ko2han/(Ko2han+So2BF(i))).*(Sno2BF(i+1)+Sno3BF(i))./(Kno3han+Sno2BF(i)+Sno3BF(i))).*fhan(i);
    SseBF(i)=DSse.*(SseBF(i+1)-2.*SseBF(i)+(SseBF(i-1)))./(hz.*i)^2 +i.*UL(i).*(SseBF(i+1)-SseBF(i))./L - (u6*nhan.*(1/YNo2han).*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno2BF(i)./(Kno2han+Sno2BF(i))).*(Ko2han./(Ko2han+So2BF(i)))).*fhan(i) -(u6*nhan.*(1/YNo3han).*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno3BF(i)./(Kno3han+Sno3BF(i))).*(Ko2han./(Ko2han+So2BF(i)))).*fhan(i) - (u6/Yhaer).*(SseBF(i)./(KSsehan+SseBF(i))).*(So2BF(i)./(Ko2han+So2BF(i))).*fhan(i);
    Sno3BF(i)=DSno3.*(Sno3BF(i+1)-2.*Sno3BF(i)+(Sno3BF(i-1)))./(hz.*i)^2 +i.*UL(i).*(Sno3BF(i+1)-Sno3BF(i))./L - (u6*nhan*((1-YNo3han)/(2.86*YNo3han)).*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno3BF(i)./(Kno3han+Sno3BF(i))).*(Ko2han./(Ko2han+So2BF(i)))).*fhan(i) + (u5/1.14).*(Snh4BF(i)./(Knh4amx+Snh4BF(i))).*(Sno2BF(i)./(Kno2amx+Sno2BF(i))).*(Ko2amx./(Ko2amx+So2BF(i))).*famx(i) + (u2/Ynb*(Sno2BF(i)./(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i)) + ((u3/Ynsp).*(Sno2BF(i)./(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i)) - ((1-fI)/2.86)*baob*naob.*(Ko2aob./(Ko2aob+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3aob+Sno2BF(i)+Sno3BF(i)).*faob(i) - ((1-fI)/2.86)*bnb*nnb.*(Ko2nb/(Ko2nb+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nb+Sno2BF(i)+Sno3BF(i)).*fnb(i) - ((1-fI)/2.86)*bnsp*nnsp.*(Ko2nsp/(Ko2nsp+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nsp+Sno2BF(i)+Sno3BF(i)).*fnsp(i) -  ((1-fI)/2.86)*bamx*namx.*(Ko2amx/(Ko2amx+So2BF(i))).*(Sno2BF(i)+Sno3(i))./(Kno3amx+Sno2BF(i)+Sno3BF(i)).*famx(i)-((1-fI)/2.86)*bhan*nhan.*(Ko2han/(Ko2han+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3han+Sno2BF(i)+Sno3BF(i)).*fhan(i);
    Sno2BF(i)=DSno2.*(Sno2BF(i+1)-2.*Sno2BF(i)+(Sno2BF(i-1)))./(hz.*i)^2 +i.*UL(i).*(Sno2BF(i+1)-Sno2BF(i))./L - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan.*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno2BF(i)./(Kno2han+Sno2BF(i))).*(Ko2han/(Ko2han+So2BF(i))).*fhan(i)) - (((u5/Yamx)+(u5/1.14)).*(Snh4BF(i)./(Knh4amx+Snh4BF(i))).*(Sno2BF(i)/(Kno2amx+Sno2BF(i))).*(Ko2amx/(Ko2amx+So2BF(i)))).*famx(i) + ((u1/Yaob).*(So2BF(i)/(Ko2aob+So2(i))).*(Snh4BF(i)./(Knh4aob+Snh4BF(i)))).*faob(i) - ((u2/Ynb).*(Sno2BF(i)./(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i)) - ((u3/Ynsp).*(Sno2BF(i)./(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i)))).*fnsp(i);
    Snh4BF(i)=DSnh4.*(Snh4BF(i+1)-2.*Snh4BF(i)+(Snh4BF(i-1)))./(hz.*i)^2 +i.*UL(i).*(Snh4BF(i+1)-Snh4BF(i))./L - (u1*(INBM+1/Yaob).*(So2BF(i)./(Ko2aob+So2BF(i))).*(Snh4BF(i)./(Knh4aob+Snh4BF(i)))).*faob(i) -(u5*(INBM+1/Yamx).*(Snh4BF(i)./(Knh4amx+Snh4BF(i))).*(Sno2BF(i)./(Kno2amx+Sno2BF(i))).*(Ko2amx./(Ko2amx+So2BF(i)))).*famx(i) - (u2*INBM.*(Sno2BF(i)./(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i)))).*(i) -(u3*INBM.*(Sno2BF(i)./(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i)))).*fnsp(i) - INBM*u6*nhan.*((SseBF(i)./(KSsehan+SseBF(i))).*(Sno2BF(i)./(Kno2han+Sno2BF(i))).*(Ko2han/(Ko2han+So2BF(i)))).*fhan(i) - u6*INBM*nhan.*((SseBF(i)./(KSsehan+SseBF(i))).*(Sno3BF(i)./(Kno3han+Sno3BF(i))).*(Ko2han./(Ko2han+So2BF(i)))).*fhan(i) - u6*INBM.*((SseBF(i)./(KSsehan+SseBF(i))).*(So2BF(i)./(Ko2han+So2BF(i)))).*fhan(i) + (INBM-(fI*INXI))*baob*naob.*((Ko2aob/(Ko2aob+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3aob+Sno2BF(i)+Sno3BF(i))).*faob(i) + (INBM-(fI*INXI))*bnb*nnb.*((Ko2nb/(Ko2nb+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nsp+Sno2BF(i)+Sno3BF(i))).*fnb(i) +(INBM-(fI*INXI))*bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nsp+Sno2BF(i)+Sno3BF(i))).*fnsp(i) + (INBM-(fI*INXI))*bamx*namx.*((Ko2amx./(Ko2amx+So2BF(i))).*(Sno2BF(i+1)+Sno3BF(i))./(Kno3amx+Sno2BF(i)+Sno3BF(i))).*famx(i) + (INBM-(fI*INXI))*bhan*nhan.*((Ko2han/(Ko2han+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3han+Sno2BF(i)+Sno3BF(i))).*fhan(i) +(INBM-(fI*INXI))*baob.*(So2BF(i)./(Ko2aob+So2BF(i))).*faob(i) + (INBM-(fI*INXI))*bnb.*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i) +(INBM-(fI*INXI))*bnsp.*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i) + (INBM-(fI*INXI))*bamx.*(So2BF(i)./(Ko2amx+So2BF(i))).*famx(i) + (INBM-(fI*INXI))*bhan.*(So2BF(i)./(Ko2han+So2BF(i))).*fhan(i);
    So2BF(i)=DSo2.*(So2BF(i+1)-2.*So2BF(i)+(So2BF(i-1)))./(hz.*i)^2 +i.*UL(i).*(So2BF(i+1)-So2BF(i))./L - -(((1-Yhaer)/Yhaer)*u6.*(SseBF(i)./(KSsehan+SseBF(i))).*(So2BF(i)./(Ko2han+So2BF(i)))).*fhan(i) -(((3.43-Yaob)/Yaob)*u1.*(So2BF(i)./(Ko2aob+So2BF(i))).*(Snh4BF(i)./(Knh4aob+Snh4BF(i)))).*faob(i) - (((1.14-Ynb)/Ynb)*u2.*(Sno2BF(i)/(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i)))).*fnb(i) - (((1.14-Ynsp)/Ynsp)*u3.*(Sno2BF(i)/(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i)))).*fnsp(i) - (1-fI)*baob.*(So2BF(i)./(Ko2aob+So2BF(i))).*faob(i) - (1-fI)*bnb.*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb -(1-fI)*bnsp.*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i) - (1-fI)*bamx.*(So2BF(i)./(Ko2amx+So2BF(i))).*famx(i) - (1-fI)*bhan.*(So2BF(i)./(Ko2han+So2BF(i))).*fhan(i);
    
    faob(i)=(muaob(i)-muaver(i)).*faob(i) - (U-UL(i).*UL(i)).*(faob(i+1)-faob(i))./L;
    fnb(i)=(munb(i)-muaver(i)).*fnb(i) - (U-UL(i).*UL(i)).*(fnb(i+1)-fnb(i))./L;
    fnsp(i)=(munsp(i)-muaver(i)).*fnsp(i) - (U-UL(i).*UL(i)).*(fnsp(i+1)-fnsp(i))./L;
    famx(i)=(muamx(i)-muaver(i)).*famx(i) - (U-UL(i).*UL(i)).*(famx(i+1)-famx(i))./L;
    fhan(i)=(muhan(i)-muaver(i)).*fhan(i) - (U-UL(i).*UL(i)).*(fhan(i+1)-fahan(i))./L;
    muaver(i)=faobi(i).*muaob(i) + fnb(i).*munb(i) + fnsp(i).*munsp(i) + famx(i).*muamx(i) + fhan(i).*muhan(i);
    UL(i)=L.*muaver(i).*i+sigma;
    VL(i)=VR-AL.*(i + LL);
end
fval=[SseBL(i);Sno3BL(i);Sno2BL(i);Snh4BL(i);So2BL(i);Xs(i);muaob(i);munb(i);munsp(i);muamx(i);muhan(i);SseBF(i);Sno3BF(i);Sno2BF(i);Snh4BF(i);So2BF(i);faob(i);fnb(i);fnsp(i);famx(i);fhan(i);muaver(i);UL(i);VL(i)];
end

Answer 1

[Cris LaPierre]

Your for loop is not executing because, after plugging in variables, your loop counter is

for i=2:1000000:0
end
i

i =

     []

so the value of I is [], resulting in fval being empty. Check the values you are using for your loop counter and correct as necessary.

[KIPROTICH KOSGEY]

Thank you!

But now after correcting i=2:100, Its giving another error:

Unable to perform assignment because the left and right sides have a different number of elements.

Error in MBBRPDE (line 107)

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

[KIPROTICH KOSGEY]

I have also tried to re-write the code as below but it is again giving an error: too many input arguments!

z1=10^-6;

N=100;

hz=N/z1;

t=0:180*24;

ic1=zeros(N,1); ic2=zeros(N,1);ic3=zeros(N,1);ic4=zeros(N,1);ic5=zeros(N,1);ic6=zeros(N,1);ic7=zeros(N,1);ic8=zeros(N,1);ic9=zeros(N,1);

ic10=zeros(N,1);ic11=zeros(N,1);ic12=zeros(N,1);ic13=zeros(N,1);ic14=zeros(N,1);ic15=zeros(N,1);ic16=zeros(N,1);ic17=zeros(N,1);ic18=zeros(N,1);

ic19=zeros(N,1);ic20=zeros(N,1);ic21=zeros(N,1);ic22=zeros(N,1);ic23=zeros(N,1);ic24=zeros(N,1);

ic=[ic1;ic2;ic3;ic4;ic5;ic6;ic7;ic8;ic9;ic10;ic11;ic12;ic13;ic14;ic15;ic16;ic17;ic18;ic19;ic20;ic21;ic22;ic23;ic24];

%solver

%options = odeset('RelTol',1e-20,'AbsTol',1e-10,'InitialStep',1e-10, 'MStateDependence','strong');

[tSol, ySol]=ode15s(@(t,y) f(t,y), t, ic,[], N);

function fval=f(t,y,N,~)

fval=zeros(length(y),1);

SseBL=zeros(N,1); Sno3BL=zeros(N,1); Sno2BL=zeros(N,1); Snh4BL=zeros(N,1); So2BL=zeros(N,1); Xs=zeros(N,1);

muaob=zeros(N,1); munb=zeros(N,1); munsp=zeros(N,1); muamx=zeros(N,1); muhan=zeros(N,1); SseBF=zeros(N,1);

Sno3BF=zeros(N,1); Sno2BF=zeros(N,1); Snh4BF=zeros(N,1); So2BF=zeros(N,1); faob=zeros(N,1); fnb=zeros(N,1);

fnsp=zeros(N,1); famx=zeros(N,1); fhan=zeros(N,1); muaver=zeros(N,1); UL=zeros(N,1); VL=zeros(N,1);

%Define the three variables

SseBL=y(1);

Sno3BL=y(2);

Sno2BL=y(3);

Snh4BL=y(4);

So2BL=y(5);

Xs=y(6);

muaob=y(7);

munb=y(8);

munsp=y(9);

muamx=y(10);

muhan=y(11);

SseBF=y(12);

Sno3BF=y(13);

Sno2BF=y(14);

Snh4BF=y(15);

So2BF=y(16);

faob=y(17);

fnb=y(18);

fnsp=y(19);

famx=y(20);

fhan=y(21);

muaver=y(22);

UL=y(23);

VL=y(24);

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;

if t>=0 && t<5

Snh4o=916; Qo=0.0452*VL*10^3/(916*24);

elseif t>=5 && t<9

Snh4o=1078; Qo=0.073*VL*10^3/(916*24);

elseif t>=9 && t<12

Snh4o=909;Qo=0.132*VL*10^3/(909*24);

elseif t>=12 && t<=15

Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);

elseif t>15 && t<=21

Snh4o=937;Qo=0.185*VL*10^3/(937*24);

elseif t>21 && t<=33

Snh4o=937;Qo=0.2286*VL*10^3/(937*24);

elseif t>33 && t<=35

Snh4o=937;Qo=0.3095*VL*10^3/(937*24);

elseif t>35 && t<=37

Snh4o=937;Qo=0.3158*VL*10^3/(937*24);

elseif t>37 && t<=49

Snh4o=937;Qo=0.348*VL*10^3/(937*24);

elseif t>49 && t<=77

Snh4o=937;Qo=0.448*VL*10^3/(937*24);

elseif t>77 && t<91

Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);

elseif t>=91 && t<=119

Snh4o=959;Qo=0.4582*VL*10^3/(959*24);

elseif t>119 && t<=127

Snh4o=959;Qo=0.1738*VL*10^3/(959*24);

elseif t>127 && t<=134

Snh4o=959;Qo=0.6246*VL*10^3/(959*24);

elseif t>134 && t<=140

Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);

elseif t>140 && t<=148

Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);

elseif t>148 && t<=155

Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);

elseif t>155 && t<=161

Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);

elseif t>161 && 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=80;

AL=120400; LL=1;

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;

%% PDEs

for i=2:N

%BC-1

SseBF(1)=0; Sno3BF(1)=0; Sno2BF(1)=0; Snh4BF(1)=0; So2BF(1)=0; faob(1)=0.4;fnb(1)=0.1;fnsp(1)=0.1;famx(1)=0.3;fhan(1)=0.1;

%BC-2

SseBF(end)=L*DLSse.*(SseBL(end)-SseBF(end))./LL; Sno3BF(end)=L*DLSno3.*(Sno3BL(end)-Sno3BF(end))/LL; Sno2BF(end)=L*DLSno2.*(Sno2BL(end)-Sno2BF(end))/LL; Snh4BF(end)=L*DLSnh4.*(Snh4BL(end)-Snh4BF(end))/LL;...

So2BF(end)=L*DLSo2.*(So2BL(end)-So2BF(end))./LL; SseBF(end)=SseBL(end); Sno3BF(end)=Sno3BL(end); Sno2BF(end)=Sno2BL(end); Snh4BF(end)=Snh4BL(end); So2BF(end)=So2BL(end);

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

Sno3BL(i)=Qo.*(Sno3o-Sno3BL(i))/VL(i)-AL.*((DSno3/LL)-UL(i)).*(Sno3BL(i)-Sno3BF(i));

Sno2BL(i)=Qo.*(Sno2o-Sno2BL(i))/VL(i)-AL.*((DSno2/LL)-UL(i)).*(Sno2BL(i)-Sno2BF(i));

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

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

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

muaob(i)=u1.*(So2BF(i)./(Ko2aob+So2BF(i))).*(Snh4BF(i)./(Knh4aob+Snh4BF(i))).*faob(i) - baob.*(So2BF(i)./(Ko2aob+So2BF(i))).*faob(i) - baob.*naob.*((Ko2aob./(Ko2aob+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3aob+Sno2BF(i)+Sno3BF(i))).*faob(i);

munb(i)=u2.*(Sno2BF(i)./(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i) - bnb.*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i) - bnb*nnb.*((Ko2nb./(Ko2nb+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nb+Sno2BF(i)+Sno3BF(i))).*fnb(i);

munsp(i)=u3.*(Sno2BF(i)./(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i) - bnsp.*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i) - bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nsp+Sno2BF(i)+Sno3BF(i))).*fnsp(i);

muamx(i)=u5.*((Snh4BF(i)./(Knh4amx+Snh4BF(i))).*(Sno2BF(i)./(Kno2amx+Sno2BF(i))).*(Ko2amx/(Ko2amx+So2BF(i)))).*famx(i) - bamx.*(So2BF(i)./(Ko2amx+So2BF(i))).*famx(i) - bamx.*namx.*((Ko2amx./(Ko2amx+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3amx+Sno2BF(i)+Sno3BF(i))).*famx(i);

muhan(i)=(u6.*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno2BF(i)./(Kno2han+Sno2BF(i))).*(Ko2han/(Ko2han+So2BF(i))) + u6.*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno3BF(i)./(Kno3han+Sno3BF(i))).*(Ko2han./(Ko2han+So2BF(i)))+u6.*((SseBF(i)./(KSsehan+SseBF(i))).*(So2BF(i)./(Ko2han+So2BF(i)))).*fhan(i))- bhan.*(So2BF(i)./(Ko2han+So2BF(i))).*fhan(i) - bhan.*nhan.*((Ko2han/(Ko2han+So2BF(i))).*(Sno2BF(i+1)+Sno3BF(i))./(Kno3han+Sno2BF(i)+Sno3BF(i))).*fhan(i);

SseBF(i)=DSse.*(SseBF(i+1)-2.*SseBF(i)+(SseBF(i-1)))./(hz.*i)^2 +i.*UL(i).*(SseBF(i+1)-SseBF(i))./L - (u6*nhan.*(1/YNo2han).*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno2BF(i)./(Kno2han+Sno2BF(i))).*(Ko2han./(Ko2han+So2BF(i)))).*fhan(i) -(u6*nhan.*(1/YNo3han).*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno3BF(i)./(Kno3han+Sno3BF(i))).*(Ko2han./(Ko2han+So2BF(i)))).*fhan(i) - (u6/Yhaer).*(SseBF(i)./(KSsehan+SseBF(i))).*(So2BF(i)./(Ko2han+So2BF(i))).*fhan(i);

Sno3BF(i)=DSno3.*(Sno3BF(i+1)-2.*Sno3BF(i)+(Sno3BF(i-1)))./(hz.*i)^2 +i.*UL(i).*(Sno3BF(i+1)-Sno3BF(i))./L - (u6*nhan*((1-YNo3han)/(2.86*YNo3han)).*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno3BF(i)./(Kno3han+Sno3BF(i))).*(Ko2han./(Ko2han+So2BF(i)))).*fhan(i) + (u5/1.14).*(Snh4BF(i)./(Knh4amx+Snh4BF(i))).*(Sno2BF(i)./(Kno2amx+Sno2BF(i))).*(Ko2amx./(Ko2amx+So2BF(i))).*famx(i) + (u2/Ynb*(Sno2BF(i)./(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i)) + ((u3/Ynsp).*(Sno2BF(i)./(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i)) - ((1-fI)/2.86)*baob*naob.*(Ko2aob./(Ko2aob+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3aob+Sno2BF(i)+Sno3BF(i)).*faob(i) - ((1-fI)/2.86)*bnb*nnb.*(Ko2nb/(Ko2nb+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nb+Sno2BF(i)+Sno3BF(i)).*fnb(i) - ((1-fI)/2.86)*bnsp*nnsp.*(Ko2nsp/(Ko2nsp+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nsp+Sno2BF(i)+Sno3BF(i)).*fnsp(i) - ((1-fI)/2.86)*bamx*namx.*(Ko2amx/(Ko2amx+So2BF(i))).*(Sno2BF(i)+Sno3(i))./(Kno3amx+Sno2BF(i)+Sno3BF(i)).*famx(i)-((1-fI)/2.86)*bhan*nhan.*(Ko2han/(Ko2han+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3han+Sno2BF(i)+Sno3BF(i)).*fhan(i);

Sno2BF(i)=DSno2.*(Sno2BF(i+1)-2.*Sno2BF(i)+(Sno2BF(i-1)))./(hz.*i)^2 +i.*UL(i).*(Sno2BF(i+1)-Sno2BF(i))./L - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan.*(SseBF(i)./(KSsehan+SseBF(i))).*(Sno2BF(i)./(Kno2han+Sno2BF(i))).*(Ko2han/(Ko2han+So2BF(i))).*fhan(i)) - (((u5/Yamx)+(u5/1.14)).*(Snh4BF(i)./(Knh4amx+Snh4BF(i))).*(Sno2BF(i)/(Kno2amx+Sno2BF(i))).*(Ko2amx/(Ko2amx+So2BF(i)))).*famx(i) + ((u1/Yaob).*(So2BF(i)/(Ko2aob+So2(i))).*(Snh4BF(i)./(Knh4aob+Snh4BF(i)))).*faob(i) - ((u2/Ynb).*(Sno2BF(i)./(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i)) - ((u3/Ynsp).*(Sno2BF(i)./(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i)))).*fnsp(i);

Snh4BF(i)=DSnh4.*(Snh4BF(i+1)-2.*Snh4BF(i)+(Snh4BF(i-1)))./(hz.*i)^2 +i.*UL(i).*(Snh4BF(i+1)-Snh4BF(i))./L - (u1*(INBM+1/Yaob).*(So2BF(i)./(Ko2aob+So2BF(i))).*(Snh4BF(i)./(Knh4aob+Snh4BF(i)))).*faob(i) -(u5*(INBM+1/Yamx).*(Snh4BF(i)./(Knh4amx+Snh4BF(i))).*(Sno2BF(i)./(Kno2amx+Sno2BF(i))).*(Ko2amx./(Ko2amx+So2BF(i)))).*famx(i) - (u2*INBM.*(Sno2BF(i)./(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i)))).*(i) -(u3*INBM.*(Sno2BF(i)./(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i)))).*fnsp(i) - INBM*u6*nhan.*((SseBF(i)./(KSsehan+SseBF(i))).*(Sno2BF(i)./(Kno2han+Sno2BF(i))).*(Ko2han/(Ko2han+So2BF(i)))).*fhan(i) - u6*INBM*nhan.*((SseBF(i)./(KSsehan+SseBF(i))).*(Sno3BF(i)./(Kno3han+Sno3BF(i))).*(Ko2han./(Ko2han+So2BF(i)))).*fhan(i) - u6*INBM.*((SseBF(i)./(KSsehan+SseBF(i))).*(So2BF(i)./(Ko2han+So2BF(i)))).*fhan(i) + (INBM-(fI*INXI))*baob*naob.*((Ko2aob/(Ko2aob+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3aob+Sno2BF(i)+Sno3BF(i))).*faob(i) + (INBM-(fI*INXI))*bnb*nnb.*((Ko2nb/(Ko2nb+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nsp+Sno2BF(i)+Sno3BF(i))).*fnb(i) +(INBM-(fI*INXI))*bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3nsp+Sno2BF(i)+Sno3BF(i))).*fnsp(i) + (INBM-(fI*INXI))*bamx*namx.*((Ko2amx./(Ko2amx+So2BF(i))).*(Sno2BF(i+1)+Sno3BF(i))./(Kno3amx+Sno2BF(i)+Sno3BF(i))).*famx(i) + (INBM-(fI*INXI))*bhan*nhan.*((Ko2han/(Ko2han+So2BF(i))).*(Sno2BF(i)+Sno3BF(i))./(Kno3han+Sno2BF(i)+Sno3BF(i))).*fhan(i) +(INBM-(fI*INXI))*baob.*(So2BF(i)./(Ko2aob+So2BF(i))).*faob(i) + (INBM-(fI*INXI))*bnb.*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb(i) +(INBM-(fI*INXI))*bnsp.*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i) + (INBM-(fI*INXI))*bamx.*(So2BF(i)./(Ko2amx+So2BF(i))).*famx(i) + (INBM-(fI*INXI))*bhan.*(So2BF(i)./(Ko2han+So2BF(i))).*fhan(i);

So2BF(i)=DSo2.*(So2BF(i+1)-2.*So2BF(i)+(So2BF(i-1)))./(hz.*i)^2 +i.*UL(i).*(So2BF(i+1)-So2BF(i))./L - -(((1-Yhaer)/Yhaer)*u6.*(SseBF(i)./(KSsehan+SseBF(i))).*(So2BF(i)./(Ko2han+So2BF(i)))).*fhan(i) -(((3.43-Yaob)/Yaob)*u1.*(So2BF(i)./(Ko2aob+So2BF(i))).*(Snh4BF(i)./(Knh4aob+Snh4BF(i)))).*faob(i) - (((1.14-Ynb)/Ynb)*u2.*(Sno2BF(i)/(Kno2nb+Sno2BF(i))).*(So2BF(i)./(Ko2nb+So2BF(i)))).*fnb(i) - (((1.14-Ynsp)/Ynsp)*u3.*(Sno2BF(i)/(Kno2nsp+Sno2BF(i))).*(So2BF(i)./(Ko2nsp+So2BF(i)))).*fnsp(i) - (1-fI)*baob.*(So2BF(i)./(Ko2aob+So2BF(i))).*faob(i) - (1-fI)*bnb.*(So2BF(i)./(Ko2nb+So2BF(i))).*fnb -(1-fI)*bnsp.*(So2BF(i)./(Ko2nsp+So2BF(i))).*fnsp(i) - (1-fI)*bamx.*(So2BF(i)./(Ko2amx+So2BF(i))).*famx(i) - (1-fI)*bhan.*(So2BF(i)./(Ko2han+So2BF(i))).*fhan(i);

faob(i)=(muaob(i)-muaver(i)).*faob(i) - (U-UL(i).*UL(i)).*(faob(i+1)-faob(i))./L;

fnb(i)=(munb(i)-muaver(i)).*fnb(i) - (U-UL(i).*UL(i)).*(fnb(i+1)-fnb(i))./L;

fnsp(i)=(munsp(i)-muaver(i)).*fnsp(i) - (U-UL(i).*UL(i)).*(fnsp(i+1)-fnsp(i))./L;

famx(i)=(muamx(i)-muaver(i)).*famx(i) - (U-UL(i).*UL(i)).*(famx(i+1)-famx(i))./L;

fhan(i)=(muhan(i)-muaver(i)).*fhan(i) - (U-UL(i).*UL(i)).*(fhan(i+1)-fahan(i))./L;

muaver(i)=faob(i).*muaob(i) + fnb(i).*munb(i) + fnsp(i).*munsp(i) + famx(i).*muamx(i) + fhan(i).*muhan(i);

UL(i)=L.*muaver(i).*i+sigma;

VL(i)=VR-AL.*(i + LL);

end

fval=[SseBL(i);Sno3BL(i);Sno2BL(i);Snh4BL(i);So2BL(i);Xs(i);muaob(i);munb(i);munsp(i);muamx(i);muhan(i);SseBF(i);Sno3BF(i);Sno2BF(i);Snh4BF(i);So2BF(i);faob(i);fnb(i);fnsp(i);famx(i);fhan(i);muaver(i);UL(i);VL(i)];

end

[Torsten]

Your code doesn't make sense. Why do you supply the solution variables SseBL, Sno3BL,... in fval ? You have to return their time derivatives dSseBLdt, dSno3BLdt,... instead.

And SseBL,Sno3BL,... are arrays of size (N,1). Thus instead of

SseBL=y(1);
Sno3BL=y(2);
...

you have to use

SseBL=y(1:N);
Sno3BL=y(N+1:2*N);
...   

Same for fval. You have to supply

[dSseBLdt;dSno3BLdt,...]

, not

[dSseBLdt(i);dSno3BLdt(i),...]

.

[KIPROTICH KOSGEY]

Thank you so much! Every time I read your answers I learn alot!

Please share the link on coding for derivates just as you have explained above for further reading.

If I have parameter whose value needs to be updated every time a solution for ODE or PDE at specific time steps is found, how do i incorporate such in the system of equations? a link on documentation would be appreciated.

[Torsten]

It would be easier to give an answer if I knew the underlying mathematical problem you are trying to solve with your code.

[KIPROTICH KOSGEY]

I would really appreciate!

I have attached an image of simplified equations. Thanks in advance.

[Torsten]

Start by discretizing and solving a simple PDE, e.g.

du/dt = d^2u/dx^2,  0 <= x <= 1, 0 <= t <= 3
u(0,t) = 0
du/dx(1,t) = u - 300
u(x,0) = 600

Compare the result with those obtained by "pdepe" to check whether your code works correctly.

Then apply what you've learned to your complex system above.

You only need to solve for fi, Si and Bi. All other variables are derived quantities that you don't need to solve for with ode15s.

[KIPROTICH KOSGEY]

Hi Torsten,

I have been working on the assignment you gave me. I have read up on pdepe and tried several examples shown on youtube and given in MATLAB documentation and it worked all along. But when I applied on the above problem you gave me I get an error 'Error using pdepe Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.' I don't know why! My code below:

m=0;

tend=0.1;

L=0.01;

x=linspace(0,L,3);

t=linspace(0,tend,3);

sol = pdepe(m,@oscpde,@oscic,@oscbc,x,t);

function [c,f,s] = oscpde(x,t,u,dudx)

c = 1;

f = dudx;

s = 0;

end

function u0 = oscic(x)

u0=600;

end

function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)

pl = 1;

ql = 0;

pr = -ur+300;

qr = 1;

end

Thanks in advance.

[KIPROTICH KOSGEY]

My bad!

I have figured out the problem. The mistake was in the value for pl: I used 1 instead of ul. Thanks.

[KIPROTICH KOSGEY]

Hi Torsten,

Sorry to bother again.

for the above problem for 'dudt=d2u/dx2' that i solved with pdepe, I have tried solving using the difference method but i am getting error: 'Index exceeds the number of array elements. Index must not exceed 1.'

code is:

nz=50;

u(1,1)=0;

u(2:nz,1)=0;

u(nz+1,1)=u(end)-300;

tSpan=[0;3];

u0=u(2:nz);

[tSol,ySol]=ode45(@(t,y) fun(t,y),tSpan,u0);

plot(tSol,ySol);

function dydt=fun(t,y)

nz=50;

z=linspace(0,0.02,1);

dz=z(1)-z(2);

u(1)=0;

u(2:nz)=0;

u(nz+1)=u-300;

dudt=zeros(nz+1,1);

for i=2:nz

dudt=(u(i+1)-u(i)+u(i-1))/dz^2;

end

dydt=dudt(2:nz);

end

I know that I have bothered you enough but please spare a minute to peruse through the above and advise. Thanks again!

[Torsten]

tstart = 0;
tend = 3;
nt = 10;
xstart = 0;
xend = 1;
nx = 50;
t = linspace(tstart,tend,nt);
x = linspace(xstart,xend,nx);

% pdepe solution
m = 0;
sol = pdepe(m,@oscpde,@oscic,@oscbc,x,t);
figure(1)
plot(x,sol(1:end,:))

% Method-of-lines solution
x = x.';
xhalf = (x(1:nx-1)+x(2:nx))/2;
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
xhalf = [xhalf;(x(end)+x(end-1))/2];
y0 = arrayfun(@(i)oscic(x(i)),1:nx);
y0(1) = 0;
tSpan = linspace(tstart,tend,nt);
[tSol,ySol]=ode45(@(t,y) fun(t,y,x,xhalf,nx),t,y0);
figure(2)
plot(x(1:end-1),ySol(1:end,:))

function [c,f,s] = oscpde(x,t,u,dudx)
  c = 1;
  f = dudx;
  s = 0;
end
function u0 = oscic(x)
  u0 = 600;
end
function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)
  pl = ul;
  ql = 0;
  pr = -ur + 300;
  qr = 1;
end
function dydt=fun(t,y,x,xhalf,nx)
  ynxplus1 = y(nx-1) + (y(nx)-300)*(x(nx+1)-x(nx-1));
  y = [y;ynxplus1];
  dydt = zeros(nx,1);
  for ix = 2:nx
    dydt(ix) = ((y(ix+1)-y(ix))/(x(ix+1)-x(ix)) - (y(ix)-y(ix-1))/(x(ix)-x(ix-1)))/(xhalf(ix)-xhalf(ix-1));
  end
end

[KIPROTICH KOSGEY]

Thank you so much Torsten, you are so kind!

I have used the central difference method and removed the 'xhalf' code and I still get the same result. However, when i introduced a second equation (du2dt=d2u/dx2+u2*(1-u1-u2). IC: u2(0,x) and BC: u2(t,0)=0 and dudx(t,l)=0...

I get the error: Index in position 2 exceeds array bounds. Index must not exceed 1. What could the problem be? Your continued assistance is highly appreciated.

code:

tstart = 0; tend = 3; nt = 10; xstart = 0; xend = 1; nx = 50;

t = linspace(tstart,tend,nt); x = linspace(xstart,xend,nx);

% Method-of-lines solution

x = x.';

% Enlarge x 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,:);

y01(1)= 0; y02(1)=0; y1=[y01 y02]; tSpan = linspace(tstart,tend,nt);

[tSol,ySol]=ode45(@(t,y) fun(t,y,x,nx),t,y1);

figure(1)

plot(x(1:end-1),ySol(1:end,:))

figure(2)

%plot(x(1:end-1),ySol(1:end,:))

function u0 = oscic(x)

u0=[600;x*(x-2)];

end

function dydt=fun(t,y,x,nx)

dydt = zeros(length(y),1);

u1=y(:,1:nx);

u2=y(:,nx+1:2*nx);

u1nxplus1 = u1(nx-1) + (u1(nx)-300)*(x(nx+1)-x(nx-1));

u2nxplus1=u2(nx-1);

u1 = [u1;u1nxplus1];

u2= [u2;u2nxplus1];

for ix = 2:nx

u1(ix) = (u1(ix+1)-2*u1(ix)+u1(ix-1))/(x(ix+1)-x(ix))^2;

u2(ix)=(u2(ix+1)-2*u2(ix)+u2(ix-1))/(x(ix+1)-x(ix))^2 + u2(ix)*(1-u1(ix)-u2(ix));

end

dydt=[u1;u2];

end

[Torsten]

I have used the central difference method and removed the 'xhalf' code and I still get the same result.

This is only the case because you use a uniform grid for x.

Although you define y1 as a row vector, it's automatically converted to a column vector by the integrator. Thus

 u1=y(:,1:nx);
 u2=y(:,nx+1:2*nx);

will throw an error.

Systems of equations that arise from the discretization of diffusion equations are stiff. Don't use ode45, but ode15s to solve these problems.

[KIPROTICH KOSGEY]

sorry but how would one then navigate through this issue? If have re-written the codes using xhalf but I still get the same error. I have also used ode15s.

[Torsten]

Did you change the assignment according to what I wrote ?

 u1=y(:,1:nx);
 u2=y(:,nx+1:2*nx);
 

And resetting u1 and u2 as

 for ix = 2:nx
   u1(ix) = (u1(ix+1)-2*u1(ix)+u1(ix-1))/(x(ix+1)-x(ix))^2;
   u2(ix)=(u2(ix+1)-2*u2(ix)+u2(ix-1))/(x(ix+1)-x(ix))^2 + u2(ix)*(1-u1(ix)-u2(ix));
 end 

makes no sense at all.

 (u1(ix+1)-2*u1(ix)+u1(ix-1))/(x(ix+1)-x(ix))^2;

and

 (u2(ix+1)-2*u2(ix)+u2(ix-1))/(x(ix+1)-x(ix))^2 + u2(ix)*(1-u1(ix)-u2(ix));

are du1/dt and du2/dt, respectively, not u1 and u2.

[KIPROTICH KOSGEY]

progress here but I am not sure if i am doind the right thing!

I have editted the code for a by changing nx to nx+1 to alter the length for IC. I have also converted u1 and u2 to column vectors. However ySOL is now alternating the columns for u1 and u2 in the responses! As a result, I get errors when trying to generate the figures 'Error using plot Vectors must be the same length.'

tstart = 0; tend = 3; nt = 10; xstart = 0; xend = 1; nx = 50;

t = linspace(tstart,tend,nt); x = linspace(xstart,xend,nx);

% Method-of-lines solution

x = x.';

xhalf = (x(1:nx-1)+x(2:nx))/2;

% Enlarge x and xhalf by a ghostcell

x = [x;x(end)+(x(end)-x(end-1))];

xhalf = [xhalf;(x(end)+x(end-1))/2];

a=arrayfun(@(i)oscic(x(i)),1:nx+1,'UniformOutput',false);

b=cell2mat(a);

y01=b(1,:); y02=b(2,:);

y01(1)= 0; y02(1)=0; y1=[y01; y02]; tSpan = linspace(tstart,tend,nt);

[tSol,ySol]=ode15s(@(t,y) fun(t,y,x,xhalf,nx),t,y1);

figure(1)

plot(x(1:end),ySol(1:end,:))

%figure(2)

%plot(x(1:end-1),ySol(1:end,:))

function u0 = oscic(x)

u0=[600;x*(x-2)];

end

function dydt=fun(t,y,x,xhalf,nx)

dydt = zeros(length(y),1);

u1=y(1:nx,:)';

u2=y(nx+1:2*nx,:)';

u1nxplus1 = u1(nx-1) + (u1(nx)-300)*(x(nx+1)-x(nx-1));

u2nxplus1=u2(nx-1);

u1 = [u1';u1nxplus1'];

u2= [u2';u2nxplus1'];

for ix = 2:nx

u1(ix) = ((y(ix+1)-y(ix))/(x(ix+1)-x(ix)) - (y(ix)-y(ix-1))/(x(ix)-x(ix-1)))/(xhalf(ix)-xhalf(ix-1));

u2(ix) = ((y(ix+1)-y(ix))/(x(ix+1)-x(ix)) - (y(ix)-y(ix-1))/(x(ix)-x(ix-1)))/(xhalf(ix)-xhalf(ix-1));

end

dydt=[u1;u2];

end

[Torsten]

As I already wrote, you don't have to redefine u1 and u2 in this loop, but you must derive du1dt and du2dt from u1 and u2.

for ix = 2:nx
    u1(ix) = ((y(ix+1)-y(ix))/(x(ix+1)-x(ix)) - (y(ix)-y(ix-1))/(x(ix)-x(ix-1)))/(xhalf(ix)-xhalf(ix-1));
    u2(ix) = ((y(ix+1)-y(ix))/(x(ix+1)-x(ix)) - (y(ix)-y(ix-1))/(x(ix)-x(ix-1)))/(xhalf(ix)-xhalf(ix-1));
end 

[Torsten]

tstart = 0;  
tend = 3;   
nt = 10;   
xstart = 0;  
xend = 1;  
nx = 50;
t = linspace(tstart,tend,nt);  
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
xhalf = (x(1:nx-1)+x(2:nx))/2;
% Enlarge x and xhalf by a ghostcell 
x = [x;x(end)+(x(end)-x(end-1))];
xhalf = [xhalf;(x(end)+x(end-1))/2];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y01 = b(1,:).';  
y02 = b(2,:).';
y01(1) = 0;  
y02(1) = 0;   
y1=[y01; y02];   
[tSol,ySol]=ode15s(@(t,y) fun(t,y,x,xhalf,nx),t,y1);
figure(1)
plot(x(1:end-1),ySol(:,1:nx))

figure(2)
plot(x(1:end-1),ySol(:,nx+1:end))

function dydt=fun(t,y,x,xhalf,nx)
  dydt = zeros(length(y),1);
  u1 = y(1:nx,:);
  u2 = y(nx+1:2*nx,:);
  u1nxplus1 = u1(nx-1) + (u1(nx)-300)*(x(nx+1)-x(nx-1));
  u2nxplus1 = u2(nx-1);

  u1 =  [u1;u1nxplus1];
  u2 =  [u2;u2nxplus1];
  du1dt = zeros(nx,1);
  du2dt = zeros(nx,1);
  for ix = 2:nx
    du1dt(ix) = (u1(ix+1)-2*u1(ix)+u1(ix-1))/(x(ix+1)-x(ix))^2;
    du2dt(ix) = (u2(ix+1)-2*u2(ix)+u2(ix-1))/(x(ix+1)-x(ix))^2 + u2(ix)*(1-u1(ix)-u2(ix));
  end 
  dydt = [du1dt;du2dt];
end
function u0 =  oscic(x)
  u0 = [600;x*(x-2)];
end

[KIPROTICH KOSGEY]

Hi Torsten

Thank you so much for the assistance you have accorded me all along. I added an ode in the above simple system of equations and everything worked fine. However, when I applied in my system of equations that I captured in the attachment previously given, it is giving an error relating to the boundary conditions: Index exceeds the number of array elements. Index must not exceed 1.

Error in MBBR>fun (line 90)

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;

the left boundary conditions are ode (dfi/dt=(mui-muaver)*fi) i=aob,nb,nsp,am,han

and also pde(dsi/dx=0) i=SseBF,Sno3BF,Sno2BF,So2BF,Snh4BF

right boundary conditions are pde(L*DLi*(SLi-Si))/LL*Di, Li=SseBL,Sno3BL,Sno2BL,So2BL,Snh4BL and i=SseBF,Sno3BF,Sno2BF,So2BF,Snh4BF

I have attached the code .

%BC-1

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

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=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;

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];