Problem with solving coupled ODE and DAE equations with mass matrix (Error using daeic12 (line 77) This DAE appears to be of index greater than 1)

Hi everyone,

I am trying to solve 6 ODE equations coupled with 1 DAE one. The ODE equations have been discritized in space domain and ode15s MATLAB solver is used to solve the equations in time domain. I have defined a diagonal mass matrix including unity elements for ODEs and zero values for DAE equation. Assuming that "n" is the No. of nodes along the space, M has been defined as follows:

```

O=ones(1,6*n);

Ze=zeros(1,n);

dia=horzcat(O,Ze);

M=diag(dia);

```

Here is the important parts of my main and function m.files:

```

.

.

.

y0(1:n)=zeros(n,1); %C0

y0(n+1:2*n)=zeros(n,1); %qc0

y0(2*n+1:3*n)=Tini*ones(n,1); %T0

y0(3*n+1:4*n)=Tini*ones(n,1); %Ts0

y0(4*n+1:5*n)=zeros(n,1); %N0

y0(5*n+1:6*n)=zeros(n,1); %qN0

y0(6*n+1:7*n)=zeros(n,1); %U0

t_int=1;

tf=3000;

nt=(tf/t_int);

tspan=(1:t_int:tf);

options=odeset('Mass',M);

[t,y]=ode15s(@adsor,tspan,y0,options);

```

and

```

function dydt=adsor(t,y)

.

.

.

dydt=zeros(7*n,1);

y(1,1)=Cin(1); %BC of C

y(n,1)=y(n-1,1); %BC of C

y(2*n+1,1)=Tg; %BC of Tg

y(3*n,1)=y(3*n-1,1); %BC of Tg

y(4*n+1,1)=Cin(2); %BC of N

y(5*n,1)=y(5*n-1,1); %BC of N

y(6*n+1,1)=Q/A; %BC of U

C(1:n,1)=y(1:n,1);

qc(1:n,1)=y(n+1:2*n,1);

T(1:n,1)=y(2*n+1:3*n,1);

Ts(1:n,1)=y(3*n+1:4*n,1);

N(1:n,1)=y(4*n+1:5*n,1);

qN(1:n,1)=y(5*n+1:6*n,1);

U(1:n,1)=y(6*n+1:7*n,1);

for i=2:n-1

dqcdt(i,1)=K1(1)*(((qmC*bC*R*Tg*C(i))./((1+(bC*R*Tg*C(i)).^ntC).^(1/ntC)))-qc(i));

dCdt(i,1)=((aEC(i).*C(i+1)+aWC(i).*C(i-1)-aPC(i).*C(i))/delz)+((((eb-1)*(1-es))/eb)*rs*dqcdt(i-1));

dqNdt(i,1)=K1(1)*(((qmN*bN*R*Tg*N(i))./((1+((bN*R*Tg*N(i))).^ntN).^(1/ntN)))-qN(i));

dNdt(i,1)=((aEC(i).*N(i+1)+aWC(i).*N(i-1)-aPC(i).*N(i))/delz)+((((eb-1)*(1-es))/eb)*rs*dqNdt(i-1));

dTdt(i,1)=((aET(i).*T(i+1)+aWT(i).*T(i-1)-(aPT(i)+((((1-eb)*as*h(i))/eb)*delz)).*T(i))/(rg*Cpg*delz))+(((((1-eb)*as*h(i))/eb).*Ts(i))/(rg*Cpg));

dTsdt(i,1)=(((h(i)*as)/(rs*Cps)).*(T(i)-Ts(i)))+((1/Cps)*(delH(1)*dqcdt(i)+delH(2)*dqNdt(i)));

dUdt(i,1)=U(i)-U(i-1)+((((1-eb)*delz)/(Ct*eb))*(dqcdt(i)+dqNdt(i)));

end

dydt(1:7*n,1)=[dCdt;0 ;dqcdt;0 ;dTdt;0 ;dTsdt;0 ;dNdt;0 ;dqNdt;0 ;dUdt;0];

```

The equation located at the last line (i.e. dUdt(i,1)=...) is the DAE equation which is recognized by the defined mass matrix (please see "M" parameter defined before). I am not sure if this method is correct to specify the DAE equations! Anyway,when I run these codes, I face the following error:

Error using daeic12 (line 77)

This DAE appears to be of index greater than 1.

I checked the Jacobian and found that it is a n*n matrix at t0. There are two columns (the first and last ones) which are all zero and I think it makes the Jacobian matrix singular (as explained in https://ch.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html).

Could you please guide me how to overcome this error and fix my code? Any help is appreciated.

Best regards