请高手帮我看看程序哪有问题

zhouwendou

function dx=jeffcott(t,X)
global w;
p=0.068;%不平衡质量偏心率
r=0.210;%轴颈半径
L=0.2016;%轴承宽度
c=0.000699;%半径间隙比
m=104860/9.8;%轴承的静载荷
n=0.0178;%润滑油的黏度
M=m*c*w*c.^2/(n*L*r.^3);%无量纲质量
G=9.8/(c*w*w);
Z1=X(1);
Z2=X(2);
Z3=X(3);
Z4=X(4);
e=sqrt(Z1.^2+Z3.^2);
ee=(Z1*Z2+Z3*Z4)/e;
Q=(Z3*Z2-Z1*Z4)/e.^2;
E=6*(1-2*Q)*(2*e*Z1/((2+e.^2)*(1-e.^2))-pi*Z3/((2+e.^2)*sqrt(1-e.^2)))+(6*ee*(pi.^2*(2+e.^2)-16)*Z1/(pi*e*sqrt(1-e.^2-4*Z3)))/((2+e.^2)*(1-e.^2));
F=6*(1-2*Q)*(pi*Z1/((2+e.^2)*sqrt(1-e.^2))-2*e*Z3/((2+e.^2)*(1-e.^2)))+(6*ee*(pi.^2*(2+e.^2)-16)*Z3/(pi*e*sqrt(1-e.^2+4*Z1)))/((2+e.^2)*(1-e.^2));
dx=zeros(4,1);
dx(1)=Z3;
dx(2)=Z4;
dx(3)=-E/M+p*sin(t);
dx(4)=-F/M+p*cos(t)+G;


clear;
w=300;
period=2*pi;
k=0;
step=2*pi/100;  
    y0=[0.1 0.1 0.1 0.1];
    k=k+1;
    tspan=[0:step:200*period];
    [t,Y]=ode45(@jeffcott,tspan,y0);
    y0=Y(end,:)
    j=1;
    for i=200:300
        tspan=[0:step:period];
        [t,Y]=ode45(@jeffcott,tspan,y0);
        YY1(k,j)=Y(end,1);  
        j=j+1;              
        y0=Y(end,:);
    end
bifdata=YY1(:,end-50:end);
plot(w,bifdata,'k.','markersize',1);


只做了w=300的情况，发现
随着时间的延长，结果发散了，不知道程序哪边出了问题，请高手帮帮忙查看一下，多谢多谢！