ODE15s-Index exceeds matrix dimensions

Hello all,

I am trying to solve 10 differential equations using ODE15s. But,

I got this error.

??? Index exceeds matrix dimensions.

Can someone help me with this code and tell me the mistake?

The mfiles are below.

function f=sofc(t,x) 

global u 

global m 

f=zeros(10,1); 

m=fsolve(@algebraeq,[5.35379707129652,0.380690723666416,0.266214812909653]); 

Rload=u; 

current=m(1); 

ethaa=m(2); 

ethac=m(3); 

dA=10^-4; 

dC=0.5*10^-4; 

dE=1.8*10^-4; 

L=0.04; 

B=0.04; 

F=96485.3; 

HA=10^-3; 

HC=10^-3; 

Taircatin=950; 

Tfuelanin=950; 

vaircatin=78.98; 

vfuelanin=2.88; 

rhocp=10^6; 

deltaHR=-241830; 

betha1=3.34*10^4; 

betha2=1.03*10^4; 

alpha=25; 

Deltaanode=5*10^-5; 

Deltacathode=5*10^-5; 

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*(Tfue-lanin^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-ethaa-ethac-rhoE*dE*current/(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)*current/(L*B)+(alpha 
+cph2solid/(2*F)*current/(L*B))*(x(10)/x(8)-x(1))+(alpha+cpo2solid/ 
(4*F)*current/(L*B))*(x(5)/x(3)-x(1))); 

f(2)=R/Deltacathode*(No2-current/(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))*(current/(2*F))); 

f(7)=R/Deltaanode*(-NH2o+(1/(L*B))*(current/(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)); 

function g=algebraeq(m) 

global u 

global x 

Rload=u; 

dE=1.8*10^-4; 

L=0.04; 

B=0.04; 

F=96485.3; 

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; 

deltaGR=-175933; 

deltaSR=-57; 

Tref=1300; 

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-m(2)-m(3)-rhoE*dE*m(1)/(L*B); 

g(1)=m(1)/(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))*m(2))-exp(- 
thetacA*F/(R*x(1))*m(2))); 

g(2)=m(1)/(L*B)-gammaC*((x(2)*x(1))/(1.013*10^5)).^0.25*exp(-EC/ 
(R*x(1)))*(exp(thetaaC*F*m(3)/(R*x(1)))-exp(-thetacC*F*m(3)/ 
(R*x(1)))); 

g(3)=v-Rload*m(1); 

and then, I have this script to integrate. 

global u 

u=5; 

options=odeset('RelTol',1e-10,'AbsTol',1e-10); 

tspan=[0 100]; 
x0=[1053.96206726841,19.7782866165184,12.0239706084605,2.40128132003310,

11462.7463172666,88.1877747305364,12.2062806809450,12.0659722222222

,
10.6185407334819,12057.5693431174,5.35379707129652,0.380690723666416,0.26621481290965¬3]; 

[t x]=ode15s(@sofc,tspan,x0,options);

Answer 1

[Andrew Newell]

@Mona, if you are writing files this long you should learn how to use the debugger (try this tutorial). Your error message should give you a file and line number for where the error occurred, e.g.,

Error in ==> algebraeq at 43

You should be able to just click on the underlined location to go to that line. Then click beside the line number for that line (you should see a red dot). Finally, run the script again. It should stop at the line number and you can look at the values of the variables by hovering your mouse over them.

[Mona]

Thanks Andrew. I will check this tutorial right now.

Answer 2

[Jan]

Instead of posting large sections of unformatted code (please read the "Markup" link on this page), the complete error message would be more helpful, because it contains the number of the line, which causes the problem. In addition post this line only.

In general "dbstop if error" helps to inspect the problems: The debugger stops, when the error occurs. Then you can observethe current variables in the Command window, Workspace Browser, etc.

[Mona]

Hi Jan

Thanks for the quick response. Where is the "Markup" link? I am not sure that I got the point. I wanna folllow this in my future posts.

[Andrew Newell]

Upper right hand corner of the box below "Add an Answer".

[Jan]

Look on top of the "Add an Answer" field. Direct link: http://www.mathworks.com/matlabcentral/answers/help/markup

[Mona]

Thanks so much.