[求助]加入for循环后为何计算出错？

[chytemp]

程序虽然有些长，我想直接贴出来比较方便吧，如下：

function N=drawmesh(m4,m5)

n=2;
R=8.314;
F=96485;
T0=298;
ra=5.5*10^8;
rc=7.0*10^8;
Ea=1.0^10^5;
Ec=1.2*10^5;
ila=2.99*10^4;
ilc=2.16*10^4;
p1=0.97;
p2=0.21;
p3=0.03;
p0=1;
Eel=8*10^4;
o=3.6*10^7;
L=20*10^(-6);
s=-55.5666;
h=-248303;
k=1/100;

i=14500;
T=1073;
%T=1173, i=17000
%T=1273, i=18500
% i can be 2000A/m^2 or 14500A/m^2 (where the power density of the fuel cell attains its maximum)

g=h-T*s;

i0a=ra*(p1/p0)*(p3/p0)*exp(-Ea/(R*T));
i0c=rc*(p2/p0)^(1/4)*exp(-Ec/(R*T));

E0=-g/(n*F);
E=E0+R*T/(n*F)*log(p1*p2^(1/2)/p3);

d1=2*n*asinh(i/(2*i0a))+2*n*asinh(i/(2*i0c))-log(1-i/ila)-log(1-i/ilc)+i*n*F*L/(o*R)*exp(Eel/(R*T));
nint=R*T*d1/(n*F);

Pc=i*E-i*R*T*d1/(n*F)-(i*n*F*k)/(R*T*d1)*(E-R*T*d1/(n*F))^2;
Nc=Pc/(-i*h/(n*F));

V=E-R*T*d1/(n*F);

m3=1.1;

x=1.2;

f=@(T2)(-T0*(exp(-(T2*x-T2-log((T*x-T2)/x/(T-T2))*x*i*m4+log((T*x-T2)/x/(T-T2))*x*i*m4*Nc+log((T*x-T2)/x/(T-T2))*x*m5*T-log((T*x-T2)/x/(T-T2))*x*m5*T0)/m3/x/(-i*m4+i*m4*Nc+m5*T-m5*T0))-1)/(-exp(-(T2*x-T2-log((T*x-T2)/x/(T-T2))*x*i*m4+log((T*x-T2)/x/(T-T2))*x*i*m4*Nc+log((T*x-T2)/x/(T-T2))*x*m5*T-log((T*x-T2)/x/(T-T2))*x*m5*T0)/m3/x/(-i*m4+i*m4*Nc+m5*T-m5*T0))+x)-(T0*(lambertw(exp((-m3*T*T2*x^2+m3*T^2*x^2+m3*T2^2*x-log((T*x-T2)/x/(T-T2))*T*T2*x^2-T*T2*x^2+log((T*x-T2)/x/(T-T2))*T^2*x^2-log((T*x-T2)/x/(T-T2))*T2^2-T*T2+2*T*T2*x-log((T*x-T2)/x/(T-T2))*T^2*x+log((T*x-T2)/x/(T-T2))*T*T2-m3*T*T2-2*m3*T*T2*x+m3*T^2*x+m3*T2^2+log((T*x-T2)/x/(T-T2))*T2^2*x)/m3/(T^2*x^2-T^2*x-T*T2*x^2+T*T2+T2^2*x-T2^2))/(x-1))*x-lambertw(exp((-m3*T*T2*x^2+m3*T^2*x^2+m3*T2^2*x-log((T*x-T2)/x/(T-T2))*T*T2*x^2-T*T2*x^2+log((T*x-T2)/x/(T-T2))*T^2*x^2-log((T*x-T2)/x/(T-T2))*T2^2-T*T2+2*T*T2*x-log((T*x-T2)/x/(T-T2))*T^2*x+log((T*x-T2)/x/(T-T2))*T*T2-m3*T*T2-2*m3*T*T2*x+m3*T^2*x+m3*T2^2+log((T*x-T2)/x/(T-T2))*T2^2*x)/m3/(T^2*x^2-T^2*x-T*T2*x^2+T*T2+T2^2*x-T2^2))/(x-1))-1)/(lambertw(exp((-m3*T*T2*x^2+m3*T^2*x^2+m3*T2^2*x-log((T*x-T2)/x/(T-T2))*T*T2*x^2-T*T2*x^2+log((T*x-T2)/x/(T-T2))*T^2*x^2-log((T*x-T2)/x/(T-T2))*T2^2-T*T2+2*T*T2*x-log((T*x-T2)/x/(T-T2))*T^2*x+log((T*x-T2)/x/(T-T2))*T*T2-m3*T*T2-2*m3*T*T2*x+m3*T^2*x+m3*T2^2+log((T*x-T2)/x/(T-T2))*T2^2*x)/m3/(T^2*x^2-T^2*x-T*T2*x^2+T*T2+T2^2*x-T2^2))/(x-1))*x-x-lambertw(exp((-m3*T*T2*x^2+m3*T^2*x^2+m3*T2^2*x-log((T*x-T2)/x/(T-T2))*T*T2*x^2-T*T2*x^2+log((T*x-T2)/x/(T-T2))*T^2*x^2-log((T*x-T2)/x/(T-T2))*T2^2-T*T2+2*T*T2*x-log((T*x-T2)/x/(T-T2))*T^2*x+log((T*x-T2)/x/(T-T2))*T*T2-m3*T*T2-2*m3*T*T2*x+m3*T^2*x+m3*T2^2+log((T*x-T2)/x/(T-T2))*T2^2*x)/m3/(T^2*x^2-T^2*x-T*T2*x^2+T*T2+T2^2*x-T2^2))/(x-1)))));
T2=fzero(f,1000)

T4=-T0*(exp(-(T2*x-T2-log((T*x-T2)/x/(T-T2))*x*i*m4+log((T*x-T2)/x/(T-T2))*x*i*m4*Nc+log((T*x-T2)/x/(T-T2))*x*m5*T-log((T*x-T2)/x/(T-T2))*x*m5*T0)/m3/x/(-i*m4+i*m4*Nc+m5*T-m5*T0))-1)/(-exp(-(T2*x-T2-log((T*x-T2)/x/(T-T2))*x*i*m4+log((T*x-T2)/x/(T-T2))*x*i*m4*Nc+log((T*x-T2)/x/(T-T2))*x*m5*T-log((T*x-T2)/x/(T-T2))*x*m5*T0)/m3/x/(-i*m4+i*m4*Nc+m5*T-m5*T0))+x)

T1=T2/x
T3=x*T4

u=(T-T1)/(T-T2);
v=(T3-T0)/(T4-T0);

Nt=1-T4*(x-1)/(m3*log(v))*((i*m4*(1-Nc)-m5*(T-T0))^(-1)-log(u)/(T2*(1-1/x)))
Pt=(-i*h/(n*F)*(1-Nc)+h*m5/(n*F*m4)*(T-T0))*Nt

N=Nc+Nt*(1-Nc-m5*(T-T0)/(m4*i))
P=(-i*h/(n*F))*N

% to plot N-i-T type ezmesh(@drawmesh,[14500,17000,1073,1273]) in the command window