在观察简化lorenz系统状态随周期力变化而变化（混沌态变化到周期态）的时候发现一个很奇怪，甚至有点诡异的问题。% 这是Sun提出的超混沌lorenz系统A
% 系统方程：
% x' = 10(y-x)
% y' = (24-4*c0*sin(u))x - xz + c0*sin(u)y
% z' = xy-8/3*z
% u' = w
% c0sin(wt)是内置驱动力
% 初值[x0 y0 z0 u0] = [0 -0.01 9 0];
clc;clear all;
w = 10; %内置频率（变量）
P = 2*pi/w;
range = [0.99337:0.00001:0.99338]; %这里是设定的周期力变化范围，以观察系统状态变化
i = 1;
x0 = [0,-0.01,9,0]; %初值
tspan = [0:2*pi/1000:1000*P];
figure
for c0 = range
% c0 = 0.99338; %内置驱动力，就是在这个点上有问题0.99336~0.99338和0.99337~0.99338这个点的状态不一致，十分诡异
% y1 = x; y2 = y; y3 = z; y4 = u
% [t,y] = ode45('LorenzA',tspan,[0,-0.01,9,0,c0]); % 这是我尝试将微分方程独立写出来的方法，但并不管用
dxdt = @(t,x)[10*(x(2)-x(1)); (24-4*c0*sin(x(4)))*x(1)-x(1)*x(3)+c0*sin(x(4))*x(2); x(1)*x(2)-8/3*x(3); w]; %微分方程函数
[t,y] = ode45(dxdt,tspan,x0);
subplot(3,4,i);
plot(t(60001:end),y(60001:end,1),'b'); %这是我观察的时域图
xlabel('t');
ylabel('X');
title(['hyperchaotic Lorenz system A /',num2str(c0)]);
% plot(y(60001:end,1),y(60001:end,3),'-k'); %这是我尝试观察的X-Z相轨图，但结果一样
% xlabel('X');
% ylabel('Z');
% title(['hyperchaotic Lorenz system A /',num2str(c0)]);
i = i+1;
end
这是我的结果：
0.99335~0.99338：这里可以看到在c0 = 0.99338时系统是处于周期状态的
0.99336~0.99339：这里可以看到在c0 = 0.99338时系统是处于混沌状态的
这就十分诡异了，明明什么条件都一样，为什么for循环会直接改变系统状态？
混沌系统状态可以跟初值、步长、系统参数有关，但两次循环中我什么都没改变啊，这简直跟鬼故事一样！！！！
尝试过将c0 = 0.99338的数据取出来看过，数据 y 的初值确实相同，大概前30000个数据都相同（不点开只能看到小数点后四位）
在40000个数据就发生了变化，那就猜测是初值的微小差异（微小到在matlab显示位数之外的差异）终于导致混沌系统在足够长的时间后呈现不同的状态，那么问题就来了，这初值明明是我自己设置的啊，真的跟鬼故事一样。
急急急！！！求解答！！！问题已经描述的很详细了，代码复制可直接运行
