DAEs (ODE 15s)

I have a system of 10 differential equations and 3 algebraic equations.

I used ODE15s to integrate. But, I am not sure why the initial value of the variables defined by the algebraic equations changes in the first calculation and it is not the same as initial condition vector x0.( x(11,x(12) & x(13)) . (mfiles dynamosofcpub2RApril7 and sofcmultipdynpub2Rapril7)

I call this function [t x]=dynamosofcpub2RApril7 in command window. Do you have any suggestions? Even if I use "if" statement in dynamosofcpub2RApril7 to apply the step change in variable u at a specific time, the sysyem does not converge.

Please find the mfiles below.

function [t x]=dynamosofcpub2RApril7
global u
u=10;
M=eye(13);
M(13,13)=0;
M(12,12)=0;
M(11,11)=0;
options=odeset('Mass',M,'MassSingular','yes','RelTol',1e-6,'AbsTol',1e-6);
tspan=[0:0.1:2000];
x0=[951.314293356128,19.9949642189062,12.0282577584507,2.40556847002325,
11427.3505058679,90.2350929738687,10.0833309517251,12.0659722222222,
10.8536790016920,11470.2069812040,0.126623256318938,0.276095749286402,0.139558206138141];
[t x]=ode15s(@sofcmultipdynpub2Rapril7,tspan,x0,options);
end

function f=sofcmultipdynpub2Rapril7(t,x)
% constants
global u
f=zeros(13,1);
% Rload=5;
% cpaircat=34.4;
% cpfuelan=31.4;
% if t<=100
% %     Rload=0.014;
% Rload=5;
% else
%     Rload=10;
% end;
Rload=u;
dA=10^-4;
dC=0.5*10^-4;
dE=1.8*10^-4;
L=0.04;
B=0.04;
F=96485.3;
naircatin=3.8*10^-2;
nfuelanin=1.39*10^-3;
HA=10^-3;
HC=10^-3;
Taircatin=950;
Tfuelanin=950;
% Taircatin=1100;
% Tfuelanin=1100;
vaircatin=78.98;
vfuelanin=2.88;
rhocp=10^6;
deltaHR=-241830;
gammaA=5.7*10^7;
gammaC=7*10^9;
thetaaC=1.4;
thetacC=0.6;
thetaaA=2;
thetacA=1;
betha1=3.34*10^4;
betha2=1.03*10^4;
EA=140000;
EC=160000;
alpha=25;
Deltaanode=5*10^-5;
Deltacathode=5*10^-5;
cpH2=30;
cpH2o=44;
cpo2=36;
epsiloncat=5*10^-1;
thoucat=3;
epsilonan=5*10^-1;
thouan=3;
deltaan=5*10^-7;
deltacat=5*10^-7;
nuH2=7.07;
nuo2=16.6;
nuN2=17.9;
nuH2o=12.7;
R=8.314;
Mo2=31.994;
MH2=2.02;
MH2o=18.02;
MN2=28.02;
naircatin=3.8*10^-2;
nfuelanin=1.39*10^-3;
yo2in=0.2;
yN2in=0.8;
vaircat=vaircatin;
vfuelan=vfuelanin;
ao2=34.57;
bo2=0.1078e-2;
do2=-784586;
ah2o=34.36;
bh2o=0.0627e-2;
ch2o=0.56012e-5;
ah2=27.67;
bh2=0.3386e-2;
an2=27.17;
bn2=0.4180e-2;
yH2in=0.9;
yH2oin=0.1;
deltaGR=-175933;
deltaSR=-57;
Tref=1300;
%physical parameters variations
Do2N2=1.013*10^-7*(x(1))^1.75*(1/Mo2+1/MN2)^0.5/((x(2)*x(1)*9.86923267*10^-6+x(5)*R*9.86923267*10^-6*(1-x(4)/x(3)))*(nuo2^(1/3)+nuN2^(1/3))^2);
Dko2=97*deltacat*abs(sqrt(x(1)/Mo2));
Do2=1/(1/Do2N2+1/Dko2);
ko2cat=(epsiloncat/(thoucat*Deltacathode))*Do2;
Ho2catin=ao2*(Taircatin-298.15)+bo2*(Taircatin^2/2-298.15^2/2)-do2*(1/Taircatin-1/298.15);
HN2catin=an2*(Taircatin-298.15)+bn2*(Taircatin^2/2-298.15^2/2);
Haircatin=yo2in*Ho2catin+yN2in*HN2catin;
Ho2cat=ao2*(x(5)/x(3)-298.15)+bo2*((x(5)/x(3))^2/2-298.15^2/2)-do2*(1/(x(5)/x(3))-1/298.15);
HN2cat=an2*(x(5)/x(3)-298.15)+bn2*((x(5)/x(3))^2/2-298.15^2/2);
Haircat=x(4)/x(3)*Ho2cat+(1-x(4)/x(3))*HN2cat;
DkH2=97*deltaan*abs(sqrt(x(1)/MH2));
DH2H2o=1.013*10^-7*x(1)^1.75*(1/MH2+1/MH2o)^0.5/((x(6)*x(1)*9.86923267*10^-6+x(7)*x(1)*9.86923267*10^-6)*(nuH2^(1/3)+nuH2o^(1/3))^2);
DH2=1/(1/DH2H2o+1/DkH2);
kH2an=(epsilonan/(thouan*Deltaanode))*DH2;
HH2anin=ah2*(Tfuelanin-298.15)+bh2*(Tfuelanin^2/2-298.15^2/2);
HH2oanin=ah2o*(Tfuelanin-298.15)+bh2o*(Tfuelanin^2/2-298.15^2/2)+ch2o*(Tfuelanin^3/3-298.15^3/3)-241814;
Hfuelanin=yH2in*HH2anin+yH2oin*HH2oanin;
HH2an=ah2*(x(10)/x(8)-298.15)+bh2*((x(10)/x(8))^2/2-298.15^2/2);
HH2oan=ah2o*(x(10)/x(8)-298.15)+bh2o*((x(10)/x(8))^2/2-298.15^2/2)+ch2o*((x(10)/x(8))^3/3-298.15^3/3)-241814;
Hfuelan=x(9)/x(8)*HH2an+(1-x(9)/x(8))*HH2oan;
DkH2o=97*deltaan*abs(sqrt(x(1)/MH2o));
DH2oH2=DH2H2o;
DH2o=1/(1/DH2oH2+1/DkH2o);
kH2oan=(epsilonan/(thouan*Deltaanode))*DH2o;

%Variables definition
No2=ko2cat*(x(4)-x(2)/(R));
NH2=kH2an*(x(9)-x(6)/(R));
NH2o=kH2oan*(x(7)/(R)-x(8)*(1-x(9)/x(8)));
U0=-1/(2*F)*(deltaGR-deltaSR*(x(1)-Tref)+R*x(1)*log((x(7)*x(1))/(x(6)*x(1)*((x(2)*x(1))/(1.013*10^5))^0.5)));
rhoE=1/betha1*exp(betha2/x(1));
% algebraic equations
v=U0-x(12)-x(13)-rhoE*dE*x(11)/(L*B);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%physical properties of the solid
cph2solid=ah2+bh2*x(1);
cpo2solid=ao2+bo2*x(1)+do2*x(1)^(-2);
cpo2cat=ao2+bo2*x(1)+do2*x(1)^(-2);
cpn2cat=an2+bn2*x(5)/x(3);
cpaircat=x(4)/x(3).*cpo2cat+(1-x(4)/x(3)).*cpn2cat;
cph2an=ah2+bh2*x(10)/x(8);
cph2oan=ah2o+bh2o*x(10)/x(8)+ch2o*(x(10)/x(8))^2;
cpfuelan=x(9)/x(8)*cph2an+(1-x(9)/x(8))*cph2oan;
% cpfuelan=31.4;

%Differential equations
f(1)=1/(rhocp*(dA+dC+dE))*((-deltaHR/(2*F)-v)*x(11)/(L*B)+(alpha+cph2solid/(2*F)*x(11)/(L*B))*(x(10)/x(8)-x(1))+(alpha+cpo2solid/(4*F)*x(11)/(L*B))*(x(5)/x(3)-x(1)));
f(2)=R/Deltacathode*(No2-x(11)/(L*B*4*F));
f(3)=1/(L*B*HC)*(naircatin-vaircat*x(3)*B*HC-No2*L*B);
f(4)=1/(L*B*HC)*(naircatin*yo2in-vaircat*x(4)*B*HC-No2*L*B);
f(5)=1/(L*B*HC*cpaircat)*(naircatin*Haircatin-vaircat*x(3)*Haircat*B*HC+L*B*alpha*(x(1)-x(5)/x(3))-No2*Ho2cat*L*B);
f(6)=R/Deltaanode*(NH2-(1/(L*B))*(x(11)/(2*F)));
f(7)=R/Deltaanode*(-NH2o+(1/(L*B))*(x(11)/(2*F)));
f(8)=1/(L*B*HA)*(nfuelanin-vfuelan*x(8)*B*HA+L*B*(NH2o-NH2));
f(9)=1/(L*B*HA)*(nfuelanin*yH2in-vfuelan*x(9)*B*HA-NH2*L*B);
f(10)=1/(L*B*HA*cpfuelan)*(nfuelanin*Hfuelanin-vfuelan*x(8)*Hfuelan*B*HA+alpha*(x(1)-x(10)/x(8))*L*B+L*B*(HH2oan*NH2o-HH2an*NH2));
f(11)=x(11)/(L*B)-gammaA*((x(6)*x(1))/(1.013*10^5))*((x(7)*x(1))/(1.013*10^5))*exp(-EA/(R*x(1)))*(exp(thetaaA*F/(R*x(1))*x(12))-exp(-thetacA*F/(R*x(1))*x(12)));
f(12)=x(11)/(L*B)-gammaC*((x(2)*x(1))/(1.013*10^5)).^0.25*exp(-EC/(R*x(1)))*(exp(thetaaC*F*x(13)/(R*x(1)))-exp(-thetacC*F*x(13)/(R*x(1))));
f(13)=v-Rload*x(11);
end

[Oleg Komarov]

format the code

[Rafael Freire]

At least format the code for readability, looking that is a kind of discouraging.