解方程结果为NaN？

ukman

有一微分方程组，,如下,x1,alph,z3为变量，m1.m2,m3,J2,k1,k2,k3,c1,c2,c3,e,tha为常量，zb,zzb,zzzb为系统激励。
m3*z3''+c3*(z3'-zb'+e*cos(tha)*alph')+k3(z3-zb+e*tha*cos(tha))=0;
m1*x1''+m2(x1''+e*sin(tha)*(alph)'')+m3*x1''+c1*x1'+k1*x1=0;
m2*e*cos(tha)*zb''+m2*e*sin(tha)*x1''+J2*(alph)''+c3*e*cos(tha)*(z3-zb+e*cos(tha)*(alph)')+k3*e*cos(tha)*(z3-zb+e*cos(tha)*(alph))+c2*(alph)'+k2*(alph)=0
用ode45解方程组，
m1=0.001;m2=21;m3=41;J2=0.18;k1=35191;k2=116;k3=40352;c1=677;c2=6;c3=495;e=0.08;tha=1.13;
ic=[0,0,0,0,0,0];
[T,Y]=ode45(@nnn_gz,[0:0.005:30],ic);

nnn_gz为

function dy=nnn_gz(t,y)
     zb = 0.01*cos(1.99*pi*t^2+2*pi*0.05*t);
     %chirp signal z0=sin(2*pi*f0*t+pi*k*t^2)
     %f0=0.05Hz;k=1.99;f1=60Hz,t=30s
     zzb=-1/100*sin(199/100*pi*t^2+1/10*pi*t)*(199/50*pi*t+1/10*pi);
     zzzb=-1/100*cos(199/100*pi*t^2+1/10*pi*t)*(199/50*pi*t+1/10*pi)^2-199/5000*sin(199/100*pi*t^2+1/10*pi*t)*pi;
     
m1=0.001;m2=21;m3=41;J2=0.18;k1=35191;k2=116;k3=40352;c1=677;c2=6;c3=495;e=0.08;tha=1.13;

B1=c3/m3;B2=k3/m3;B3=B1;B4=B2;A1=c3*e*cos(tha);A2=k3*e*cos(tha);B5=A1/m3;B6=A2/m3;
A3=m1+m2+m3;A4=m2*e*sin(tha);A5=c2+c3*e^2*(cos(tha))^2;A6=k2+k3*e^2*(cos(tha))^2;
A7=k3*e*cos(tha)+c3*e*cos(tha);A8=m2*e*cos(tha);A9=A7;
D1=-A4/A3;D2=-c1/A3;D3=-k1/A3;
E1=-A8/(A4*D1+J2);E2=A9/(A4*D1+J2);E3=A4*D2/(A4*D1+J2);E4=-A4*D3/(A4*D1+J2);
E5=-A6/(A4*D1+J2);E6=-A5/(A4*D1+J2);E7=-A7/(A4*D1+J2);
F1=D1*E3+D2;F2=D1*E4+D3;F3=D1*E1;F4=D1*E2;F5=D1*E5;F6=D1*E6;F7=D1*E7;
%y(1)=x1;y(2)=dx1;y(3)=alph,y(4)=dalph;y(5)=z3;y(6)=dz3;
dy=zeros(6,1);
dy(1)=y(2);
dy(2)=F1*y(2)+F2*y(1)+F3*zzzb+F4*zb+F5*y(3)+F6*y(4)+F7*y(5);
dy(3)=y(4);
dy(4)=E1*zzzb+E2*zb+E3*y(2)+E4*y(1)+E5*y(3)+E6*y(4)+E7*y(5);
dy(5)=y(6);
dy(6)=B1*zzb+B2*zb+B3*y(6)+B4*y(5)+B5*y(4)+B6*y(3);
end
计算结果Y(:,1)为NaN，数量级达到10e302，请大家帮忙看看，问题很长，不胜感激！！！

ukman

已经查出来了，是降阶的时候出了问题