如何更改为求刚度比的分岔图程序？

Rotor2014

如何更改为求解刚度比的分岔图，关键是看不懂如果改为刚度比的话，步长怎么处理？
clear all;
clc;
global w;
global u;
%global E;
%global beta;
u=0.00007;
f0=250;
range=[2.5:0.005:4];
k=1;
YY=[];%初始化YY
z0=[0 0.001 1.27e-3 0.001];%设置初值
for w=range
period=2*pi/w;
step=period/512;
w
%j=j+1;
% discard the first *** periodic data;
%除去前面2000个周期的数据，并将最后的结果作为下一次积分的初值
tspan=[0:step:2000*period];
%options=odeset('RelTol',10^-3,'AbsTol',10^-5);
[t,Y]=ode45('rubbing',tspan,z0);
z0=Y(end,:)
%j=1;
%for i=2000:2200
  tspan=[2000*period:step:2100*period];
  %options=odeset('RelTol',10^-3,'AbsTol',10^-5);
  [t,Y]=ode45('rubbing',tspan,z0);
  YY(k,:)=Y(1:512:end,1); % getthe omega data from every period end
  %j=j+1;   
  z0=Y(end,:);
  k=k+1;
%end
end
bifdata=YY(:,end-51:end);
plot(range,bifdata,'k.','LineWidth',1);
xlabel('频率比w/w0');ylabel('x');title('随频率比变化的分岔图');
%xlabel('刚度比kc/k0');ylabel('x');title('随刚度比变化的分岔图');