转子动力学相图的问题，程序调试不出来，希望高手指点

wangboy

function wwww

fs=1000;
ff=2*pi/fs;

tf=1000;
tspan=0:ff:tf;
xo=[0.001,0.001,0.001,0.01,0.001,0.001,0.001,0.001,0.001,0.01];
options=odeset('RelTol',1e-8,'AbsTol',1e-8);
[t,y]=ode15s(@w,tspan,xo,options);

plot(y(end-40000:end-5000,1),y(end-40000:end-5000,3))


function dx=w(t,x)
R=0.215;
m1=7;
m2=2;
s=2.2;
l=0.345;Re=0.2;Rd=0.12;
r=0.0175;R1=Re+Rd;d=2*r;
delta2=0.6e-3;
e1=0.06e-3;
d1=0.2;
d2=0.4*sqrt(m1/m2);
J=m1*R^2/2;I=pi*d^4/64;
E=2e11;
k=12*E*I/l^3;
w=s*sqrt(k/m1);
cT=2;
alf=26.6;belt=30;
ga=0;  
v=(w)*R1;         
F=308.53-4.13*alf-3.24*belt+0.091*alf^2+0.047*belt^2+0.043*ga^2+0.107*v^2;
belt1=belt*pi/180;
T=2*pi/w;
mm=fix(t/(T/4));
tt=t-mm*T/4;
syms theta
RR=0.02; L=0.008; miu=0.01;  dert=0.00011;

eb=sqrt(x(5)^2+x(7)^2);
cs=x(5)/eb;
ss=x(7)/eb;
deb=x(6)*cs+x(8)*ss;
dphab=(-x(6)*ss+x(8)*cs)/eb;

f1=(cos(theta))^2/(1+eb*cos(theta))^3;
f2=cos(theta)*sin(theta)^2/(1+eb*cos(theta))^3;
f3=(sin(theta))^2/(1+eb*cos(theta))^3;
theta1=atan(-deb/eb/dphab);

I1=int(f1,'theta',theta1,theta1+pi);
I2=int(f2,'theta',theta1,theta1+pi);
I3=int(f3,'theta',theta1,theta1+pi);

Fr=deb*I1+dphab*I2*eb^2;
Ft=deb*I2+dphab*I3*eb^2;
cc=miu*RR*L^3/dert^2;
Fx=cc*Fr*cs-cc*Ft*ss;
Fy=cc*Fr*ss+cc*Ft*cs;


dx=zeros(10,1);
x(5)
dx(1)=x(2);
dx(2)=F*abs(sin(2*t))*sin(tt+belt1)/m1/delta2/w^2-1/s^2*(x(1)-x(5))-2*d1/s*x(2)+e1/delta2*cos(t)+pi/2*dx(10)*e1*w^2/delta2*sin(t);
dx(3)=x(4);
dx(4)=F*abs(sin(2*t))*cos(tt+belt1)/m1/delta2/w^2-1/s^2*(x(3)-x(7))-2*d1/s*x(4)+e1/delta2*cos(t)-pi/2*dx(10)*e1*w^2/delta2*sin(t);
dx(5)=x(6);
dx(6)=Fx/m2/delta2/w^2-1/s^2*(x(1)-x(5))-2*d2/s*x(6);
dx(7)=x(8);
dx(8)=Fy/m2/delta2/w^2-1/s^2*(x(3)-x(7))-2*d2/s*x(8);
dx(9)=x(10);
dx(10)=F*abs(sin(2*t))*R1/(J*w^2*pi/2)-cT/J/w*x(10)+m1*e1*dx(2)*delta2/(J*pi/2)*sin(t)-m1*e1*dx(4)*delta2/(J*pi/2)*cos(t)-2*cT/J/w/pi;