Unable to get the following DAE system of equations solved

%% DAE of solution
tSpan = [0.000001,1];
Y0= [0.02646;0.97354;1;101325000;0.0001;0.0001;0.0001;0.0001];
MassM=[1 0 0 0 0 0 0 0;
    0 1 0 0 0 0 0 0;
    0 0 1 0 0 0 0 0;
    0 0 0 1 0 0 0 0;
    0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0];
opt=odeset('Mass', MassM);
[tSol,YSol]=ode15i(@(t,Y) YTHFcountercurrentflowFun(t,Y), tSpan, Y0, opt);

Error using solution>@(t,Y)YTHFcountercurrentflowFun(t,Y)
Too many input arguments.

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15i (line 116)
odearguments(FcnHandlesUsed, false, solver_name, ode, tspan, y0, ...

%% Plot solutions
plot(tSol, YSol);
function dY = YTHFcountercurrentflowFun (t, Y)
% CROSSFLOW cross flow model
global c1 c2 z ph pl permN2 d do Nf Uf mu1 R T di pi
c1=0.00000000004;
c2=0.0004696;
d=0.0023876;
ph= 101325;
pl= 0.001;
permN2= 0.00000000000000104836;
do= 0.0002; 
z = 0.2;
Nf= 10000;
Uf= 0.1;
mu1= 0.00001805;
R= 8.314;
T= 303.15;
di= 0.0001;
pi= 3.14159265;
s= (pi*do*z*Nf*permN2*ph)/(d*Uf);
A= (256*mu1*R*T*Uf*Uf*d)/(pi*pi*permN2*ph*ph*ph*do*di*di*di*di*Nf*Nf);
% Functions
dY(1,1)= (1/Y(3))*(Y(8)/permN2)*(Y(1)-(sqrt(Y(4))*Y(5)))*((Y(1)/Y(5))-1);
dY(2,1)= (1/Y(3))*(Y(2)-(sqrt(Y(4))*Y(6)))*((Y(2)/Y(6))-1);
dY(3,1)= (-(1-Y(4))/Y(7));
dY(4,1)= A* (1- Y(3));
dY(5,1)= Y(5)-(((Y(8)/permN2))*Y(1)*Y(7))/(1-sqrt(Y(4))+(Y(7)*(Y(8)/permN2)*sqrt(Y(4))));
dY(6,1)= Y(6)-(Y(2)*Y(7))/(1-sqrt(Y(4))+(Y(7)*sqrt(Y(4))));
dY(7,1)= Y(7)- Y(1)+ ((Y(8)/permN2)*Y(2));
dY(8,1)= Y(8)-c1*(exp(c2*Y(1)*(sqrt(Y(4))/pl)));
end

Answer 1

[Walter Roberson]

ode15i passes three parameters to the given function: (t,y,y′)

You are not required to do anything with the third parameter... but if you do not do anything with the third parameter then it is a waste to use ode15i.

Answer 2

[Torsten]

The time derivative of the first equation is 1e45 at the beginning. Check your parameters and their units.

%% DAE of solution
tSpan = [0.000001,1];
Y0= [0.02646;0.97354;1;101325000];
[tSol,YSol]=ode15s(@(t,Y) YTHFcountercurrentflowFun(t,Y), tSpan, Y0)
%% Plot solutions
plot(tSol, YSol(:,2));

function dY = YTHFcountercurrentflowFun (t, Y)
% CROSSFLOW cross flow model
c1=0.00000000004;
c2=0.0004696;
d=0.0023876;
ph= 101325;
pl= 0.001;
permN2= 0.00000000000000104836;
do= 0.0002; 
z = 0.2;
Nf= 10000;
Uf= 0.1;
mu1= 0.00001805;
R= 8.314;
T= 303.15;
di= 0.0001;
pi= 3.14159265;
s= (pi*do*z*Nf*permN2*ph)/(d*Uf);
A= (256*mu1*R*T*Uf*Uf*d)/(pi*pi*permN2*ph*ph*ph*do*di*di*di*di*Nf*Nf);
Y(8) = -(-c1*(exp(c2*Y(1)*(sqrt(Y(4))/pl))));
Y(7) = -(- Y(1)+ ((Y(8)/permN2)*Y(2)));
Y(6) = -(-(Y(2)*Y(7))/(1-sqrt(Y(4))+(Y(7)*sqrt(Y(4)))));
Y(5) = -(-(((Y(8)/permN2))*Y(1)*Y(7))/(1-sqrt(Y(4))+(Y(7)*(Y(8)/permN2)*sqrt(Y(4)))));
% Functions
dY(1,1)= (1/Y(3))*(Y(8)/permN2)*(Y(1)-(sqrt(Y(4))*Y(5)))*((Y(1)/Y(5))-1);
dY(2,1)= (1/Y(3))*(Y(2)-(sqrt(Y(4))*Y(6)))*((Y(2)/Y(6))-1);
dY(3,1)= (-(1-Y(4))/Y(7));
dY(4,1)= A* (1- Y(3));
%dY
end