matlab 计算结果为nan,matlab 计算 结果总是为Nan

本人刚刚接触matlab，对这些运算不是很懂  计算ni的位置(exp(-E_g./(2.*k.*T))+eps)这个值之前一直是0 加了eps后就有结果显示了 后面部分exp((alfa.*T)./(2.*k.*(T+beita)))一直为无穷大然后结果就为Nan 不知道怎么解决 求大神帮忙

q = 1.602*10^(-19);

k = 1.38*10^(-23);

fai_msf = -0.08;

fai_msb = 0.17;

kexi_Si = 1.05*10^(-12);

kexi_ox = 3.45*10^(-13);

d_fox = 10^(-6);

d_box = 3.8*10^(-5);

d_Si = 3*10^(-5);

%ni = 1.45*10^10;

Dn = 34;

%ni = 1.45*10^10;

%fai_f = 0.29;

fai_T = 0.026;

N_a = 10^15;

N_d = 10^20;

tao = 10^(-7);

s_0 = 10^4;

s_dSi = 10^4;

y = 10^(-7);

Q_oxf = 3*10^10*q;

Q_oxb = 10^10*q;

Q_dep = -q*N_a*d_Si;

Q_sb = 0;

V_BG = 0;

V_A = -1;

V_GK = 0.4;

E_g = 1.21;

alfa = 4.73*10^(-4);

beita = 636;

T = 250:10:350;

ni=3.39.*10^15.*T.^(3/2).*(exp(-E_g./(2.*k.*T))+eps).*exp((alfa.*T)./(2.*k.*(T+beita)));

fai_f = (k.*T./q).*log(N_a./ni);

fai_sb = fai_f;

C_Si = kexi_Si./d_Si;

C_ox = kexi_ox./d_fox;

C_box = kexi_ox./d_box;

A = -fai_msf+Q_oxf/C_ox+Q_dep/(2*C_ox)-(C_Si/C_ox)*(C_box./(C_box+C_Si))*(fai_msb-Q_oxb./C_box-Q_dep./(2*C_box));

B = 1+C_Si./C_ox-(C_Si./C_ox)*(C_box./(C_box+C_Si))*(C_Si./C_box);

C = V_BG-fai_msb+Q_oxb./C_box;

D = C_Si/C_box;

E = Q_dep/(2*C_box);

F = 1+C_Si/C_box;

G = q*N_a*y^2/(2*kexi_Si);

H = (q*N_a*d_Si)/(2*kexi_Si);

fai_sf = (V_GK+A)./B;

fai_sb = (C+D*fai_sf+E)/F;

fai_y = G+((fai_sf-fai_sb)./d_Si-H).*y+fai_sf;

% fai_sf1 = (V_GK2+A)./B;

% fai_sb1 = (C+D*fai_sf1+E)/F;

% fai_y1 = G+((fai_sf1-fai_sb1)./d_Si-H).*y+fai_sf1;

% n_VGK2 = ni*exp((fai_y1-fai_f)/fai_T);

n_y = ni.*exp((fai_y-fai_f)./fai_T);

p_y = ni./exp((fai_y-fai_f)/fai_T);

V_D = fai_T.*log(p_y.*N_d./(ni.^2));

V_D1 = fai_T.*log(N_d.*n_y/(ni.^2));

V_D0 = fai_T.*log(N_a.*N_d./(ni.^2));

V_D10 = fai_T*log(N_d/N_a);

fai_N = (V_D0+V_D10)/2-V_A;

V_A1 = -(fai_N-(fai_y-fai_f)-V_D);

Ld = (2.*kexi_Si.*(V_D-V_A1).*(N_a+N_d)./(q.*N_a.*N_d)).^(0.5);

xp = (2.*kexi_Si.*(V_D-V_A1).*N_d./(q.*N_a.*(N_a+N_d))).^(0.5);

L_d1 = (2.*kexi_Si.*(V_D+(1+V_A1)).*(N_a+N_d)/(q.*N_a.*N_d)).^(0.5);

xp1 = (2.*kexi_Si.*(V_D+(1-V_A1)).*N_d/(q.*N_a.*(N_a+N_d))).^(0.5);

tao_eff = 1/(p_y/(tao*ni)+(s_0+s_dSi)/d_Si);

% tao_eff = tao*ni/p_y;

Ln = (Dn*tao_eff)^(0.5);

L = 2*10^(-3);

L1 = L-xp-xp1;

x = 0:L1/100:L1;

nx = (N_d*exp(q*V_A1/(k*T))-n_y)*sinh(x/Ln)/sinh(L1/Ln)+n_y;

Jnx = q*Dn*(N_d*exp(q*V_A1/(k*T))-n_y)*cosh(x/Ln)/(sinh(L1/Ln)*Ln);

figure(1),plot(x,nx);

figure(2),plot(x,Jnx);