解一阶多元微分方程问题

leohust

一组一阶9元一次微分方程组，在别人的帮助下，写出了m文件及命令文件，但是调试时总是算不出来，得到的总是y =
      NaN +    NaNi      NaN +    NaNi      NaN +    NaNi      NaN +    NaNi 不知道是什么原因，还请大家帮忙看看：
function f=odefun1(t,y,a)
%Vp
Vp=pi*y(3)^3/6;
%qpr
mc=pi*y(3)^2*y(4)*y(6)*a(1)*exp(-a(4)/y(1));
qpr=a(2)*mc;
%qrad
qrad=pi*y(3)^2*a(5)*(a(6)^4-y(1)^4);
%qconv
qconv=2*a(7)*pi*y(3)*(y(1)-y(2));
%N
Vc=pi*a(16)^3/6;
n0=6*a(15)*y(5)/(pi*y(4)*(y(3)^3+a(15)*a(16)^3));
N=n0*Vc;
%qgr
mgrv1=a(8)*y(7)^(-0.3)*y(6)^1.3*exp(-a(11)/(a(12)*y(2)));
mgrv2=a(13)*y(8)^(-0.1)*y(6)^1.85*exp(-a(11)/(a(12)*y(2)));
qgr=y(5)*Vc*(mgrv1*a(10)+mgrv2*a(14));
%qp
mv=0.3*a(17)*y(9)*exp(-a(19)/y(1))+0.6*a(18)*y(9)*exp(-a(20)/y(1));
qp=a(3)*(mv+mc)*(y(1)-a(6));
%qc
lambda=a(7);
qc=2*lambda*pi*a(16)*(a(6)-y(2));
f=[(qpr+qrad-qconv)/(y(4)*Vp*a(3));
(qgr+N*qconv+N*qp+qc)/(y(5)*Vc*a(9));
-2*mc/(pi*y(4)*y(3)^2);
-6*mv/(pi*y(3)^3);
n0*(mv+mc);
-8*n0*mc/(3*y(5))-4*mgrv1-3*mgrv2;
0.3*n0*a(17)*y(9)*exp(-a(19)/y(1))/y(5)+mgrv1;
0.6*n0*a(18)*y(9)*exp(-a(20)/y(1))/y(5)+mgrv2;
a(21)-a(17)*y(9)*exp(-a(19)/y(1))-a(18)*y(9)*exp(-a(20)/y(1))];
a=[7630 32805 1.1 17615 4.7628*10^(-11) 1273 0.05 1.5*10^7 1 3.95*10^4 125000 8.314 1.3*10^9 3.21*10^4 1 0.01 3.7*10^5
1.46*10^13 8899 30186 2.7*10^(-10)];
tspan=[0 1];
y0=[300 1273 7.5*10^(-5) 1250 1.183 0.23 0 0 2.7*10^(-10)];
options=odeset('RelTol',1e-6,'AbsTol',[1e-6],'OutputFcn',@odeplot,'OutputSel',[1]);
[t,y]=ode45(@odefun1,tspan,y0,options,a);