做过转子碰摩的前辈请进来帮忙看下程序

wangu

%非线性刚度转子系统碰摩的混沌行为
%罗跃纲.2002.东北大学学报.vol.23，no.9
function dx=rubbing(t,x);
global gamma;
dx=zeros(4,1);
%parameter
xi=0.12;
alpha=0.5;
beta=3;
f=0.12;
delta=0.16;
b=0.136;
g=9.8;
f0=25;
w0=2*pi*f0;
e=sqrt(x(1)^2+x(3)^2);
H=1/2*(sign(abs(e-1))+sign(e-1));%判断是否发生碰磨
%dimensionless equations
dx=[x(2);
    -2*xi*x(2)-x(1)-alpha*(x(1)^2+x(3)^2)*x(1)-beta*(1-1/e)*(x(1)-f*x(3))*H+b*gamma^2*cos(gamma*t);
    x(4);
    -2*xi*x(4)-x(3)-alpha*(x(1)^2+x(3)^2)*x(3)-beta*(1-1/e)*(f*x(1)+x(3))*H+b*gamma^2*sin(gamma*t)-g/(delta*w0^2)];


主程序
clear all;
clc;
global gamma;
range=[0.10:0.01:1.90];
period=2*pi;
k=0;
YY1=[];
step=period/512;  
for gamma=range
    x0=[0.05 0.005 0.05 0.005];
    gamma
    k=k+1;
     % discard the first *** periodic data;
     tspan=[0:step:2000*period];
     options=odeset('RelTol',10^-3,'AbsTol',10^-5);
    [t,x]=ode45('rubbing',tspan,x0,options);
    x0=x(end,:);
    j=1;
    for i=2000:2200
        tspan=[i*period:step:(i+1)*period];
        options=odeset('RelTol',10^-3,'AbsTol',10^-5);
        [t,x]=ode45('rubbing',tspan,x0,options);
        YY1(k,j)=x(end,1);   % get the omega data from every period end
        j=j+1;               
        x0=x(end,:);
    end
end
bifdata=YY1(:,end-51:end);
plot(range,bifdata,'k.','markersize',1);
xlabel('频率比w/w0');ylabel('x');title('随频率比变化的分岔图');

我得到的图是几根线，跟原文中的分岔图差的太多了，不知道错在了哪里？请帮忙看下

wangu

原文分岔图

dollfish000

学习了，正想涉及非线性