不解：Chen's吸引子的程序问题

xiaocheng_2007

下面是采用四阶龙格库塔法求吸引子的程序。采用的是吕金虎著《混沌时间序列分析及其应用》第18页的参数设置，可是我却得不出书中的吸引子图。请大家多多指教！
%４阶龙格库塔法求chen's吸引子
% 方程表达式
% dx/dt = a*(y-x)
% dy/dt = (c-a)*x - x*z + c*y
% dz/dt = x*y - b*z
clc;
close all;
clear all;
%参数值
a = 35;
b = 3;
c = 28;
%初始值
x_0 = 0;
y_0 = 1;
z_0 = 0;
h = 0.0001;             % 积分时间步长
step1 = 30000;        % 前面的迭代点数
step2 = 5000;         % 后面的迭代点数
X = [];
Y = [];
Z = [];
for(i = 1:1:(step1 + step2))
    %e1
    x_e1 = -a*x_0 + a*y_0;
    y_e1 = c*y_0 + (c - a)*x_0 - x_0*z_0;
    z_e1 = -b*z_0 + x_0*y_0;
    %e2
    y_h = y_0 + 0.5*h*y_e1;
    x_e2 = -a*(x_0 + 0.5*h*x_e1) + a*y_h;
   
    x_h = x_0 + 0.5*h*x_e1;
    z_h = z_0 + 0.5*h*z_e1;
    y_e2 = c*(y_0 + 0.5*h*y_e1) + (c - a)*x_h - x_h*z_h;
   
    z_e2 = -b*(z_0 + 0.5*h*z_e1) + x_h*y_h;
    %e3
    y_h = y_0 + 0.5*h*y_e2;
    x_e3 = -a*(x_0 + 0.5*h*x_e2) + a*y_h;
   
    x_h = x_0 + 0.5*h*x_e2;
    z_h = z_0 + 0.5*h*z_e2;
    y_e3 = c*(y_0 + 0.5*h*y_e2) + (c - a)*x_h - x_h*z_h;
   
    z_e3 = -b*(z_0 + 0.5*h*z_e2) + x_h*y_h;
    %e4
    y_h = y_0 + h*y_e3;
    x_e4 = -a*(x_0 + h*x_e3) + a*y_h;
   
    x_h = x_0 + h*x_e3;
    z_h = z_0 + h*z_e3;
    y_e4 = c*(y_0 + h*y_e2) + (c - a)*x_h - x_h*z_h;
   
    z_e4 = -b*(z_0 + h*z_e2) + x_h*y_h;
    %叠代
    x_1 = x_0 + 1/6*h*(x_e1 + 2*x_e2 +2*x_e3 + x_e4);
    y_1 = y_0 + 1/6*h*(y_e1 + 2*y_e2 +2*y_e3 + y_e4);
    z_1 = z_0 + 1/6*h*(z_e1 + 2*z_e2 +2*z_e3 + z_e4);
   
    X = [X,x_1];
    Y = [Y,y_1];
    Z = [Z,z_1];
   
    x_0 = x_1;
    y_0 = y_1;
    z_0 = z_1;
end
X = X(10000+1:end);
Y = Y(10000+1:end);
Z = Z(10000+1:end);
figure(1);
plot3(Z,Y,X);
grid on
xlabel('x');
ylabel('y');
zlabel('z');

yuling

我用matlab计算了一下，和书上的差不多，程序如下：
function dy = Chens(a,y,c)
a=35;
c=28;
b=3;
dy=zeros(3,1);
dy(1)=-a*(y(1)-y(2));
dy(2)=-y(1)*y(3)+(c-a)*y(1)+c*y(2);
dy(3)=y(1)*y(2)-b*y(3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
clear all
t0=0;
tf=35;
[t,x]=ode45(@Chens,[t0:0.0001:tf],[0,1,0]);
plot3(x(10000:end,3),x(10000:end,2),x(10000:end,1));

xiaocheng_2007

马虎了，已经解决。迭代数太少了。。。:@Q