如何根据如下方程（还有我的程序）画出图片中的分岔图，poincare映射图？恳请指正！

chunshui2003

方程及分岔图、映射图以图片的形式上传。

我的程序：

function y=fangcheng(t,x)
% 求解定值问题

m=20;
%转子质量
e=0.0009;
%圆盘偏心距
xi=0.13;
xi_z=0.15;
omega_0=16;
omega_e=22;
%轴向推力的频率
delta=0.0018;
%静止时转子与定子之间的间隙
kc=25*10^5;
%定子的径向刚度
f=0.1;
%转子与定子之间的摩擦系数
R=0.1;
%碰摩点到转子形心的距离
Fe=9;
%轴向推力的幅值
omega_z=20;
g=9.80;
%重力加速度

global omega


Vr=R*omega+(x(2).*x(5)-x(3).*x(4))/sqrt(x(2).^2+x(4).^2);
x(1)=t;

y=zeros(7,1);
y(1)=t;
y(2)=x(3);
y(3)=e*omega.^2*cos(2*pi*omega.*x(1))-2.*xi*omega_0.*x(3)-omega_0^2.*x(2)-kc/m*(1-delta/sqrt(x(2).^2+x(4).^2)).*(x(2)-f*Vr/sqrt(Vr^2+x(7).^2).*x(4));
y(4)=x(5);
y(5)=e*omega.^2*sin(2*pi*omega.*x(1))-2.*xi*omega_0.*x(5)-omega_0^2.*x(4)-kc/m*(1-delta/sqrt(x(2).^2+x(4).^2)).*(x(4)+f*Vr/sqrt(Vr^2+x(7).^2).*x(2));
y(6)=x(7);
y(7)=Fe/m*cos(2*pi*omega_e*x(1))-g-2*xi_z*omega_z*x(7)-omega_z^2*x(6)-kc/m*f*(sqrt(x(2)^2+x(4)^2)-delta)-x(7)/sqrt(Vr^2+x(7)^2);






%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear

global omega

tic

omega=2.51*16;
period=1/omega;
tspan=[0:period/100:200*period];
y0=[1 0.1 0.1 0.1 0.1 0.1 0.1];
[t,x]=ode45('fangcheng',tspan,y0);

toc

然后就是画分岔图和映射图了。按照2.51倍的omega_0计算，希望能够得到和图中一样的结果。但是映射图和轴心轨迹图怎么也对不上，1.39倍的也对不上，我想我是哪里错了。请问是我的计算程序有问题还是计算后画图取点有问题。

轴心轨迹图：
Plot(x(4000:end,2),x(4000:end,4));
(其他的取点也试过了，但得不到满意的结果)

映射图：
Plot(x(4000:100:end,2),x(4000:100:end,3));
（依然不行）

希望各位如果有时间帮忙看一下，到底是哪里的问题，我需要改进哪些部分。以后如果画图的时候应该注意什么，因为这些真的是属于自己琢磨，没有师兄指点，有时候很头疼，就像一道障碍很难越过。先提前谢谢大家了！

另外，分岔图的程序（试了也不行），也请大家看一下
clc
clear
% global omega
y0=[1 0.1 0 0.1 0 0.1 0];
omega=16:100;
for i=1:length(omega)

disp(omega(i))

period=1/omega(i);

tspan=[0:period/100:100*period];

[t,x]=ode45('fangcheng',tspan,y0,[],omega(i));

y0=x(end,:);

plot(omega(i),x(2000:100:end,2),'markersize',2);

hold on;
end

这里面把前面的方程改动一下，将omega作为参数输入到函数中，其余的没变。

htwei

我看不出错误  你看看是不是可以把步长减小点 period=1/omega(i);

tspan=[0:period/1000:100*period];