怎样利用４阶龙格库塔法画带脉冲的时滞系统

lyapunov

想问一下怎样利用４阶龙格库塔法画带脉冲的时滞系统
我只有画不带时滞的脉冲系统相图的程序，但是改不成加时滞的，看到一位叫何迎生的老师的一篇文章说可以利用４阶龙格库塔法画带脉冲的时滞系统的相图。纠结的是这位何老师并没有写联系方式。
                                                                                                                   小女先在此谢过各位大侠的指点


%四阶龙格库塔法求 Rank one 奇怪吸引子.
%方程表达式
%dx/dt=2*b*x(t)*y(t)+c*x(t)+c*f*x(t)*x(t)-d*x(t)
%dy/dt=1-2*a*x(t-tau)-b*y(t)*y(t)
clc;
close all;
clear all;
%参数值
a = 1;
c = -1;
b =1;
d=0.8;
f=-1.23;
T=8;
p=0.5;
%初始值
x_0 =0.07;
y_0 =0.5;
z_0=0.8;
X = [];
Y = [];
Z = [];
pn=50;               % 周期数
h =0.01;            % 积分时间步长
step1 =T/h;         % 前面的迭代点数
m=pn*T/h;
for j = 1:m-1
    k=fix(j/step1);
   
    X = [X,x_0];
    Y = [Y,y_0];
    Z = [Z,z_0];
   
       tf=(j>=k*step1&j<k*step1+p/h);
    %e1
    x_e1 =2*b*x_0*y_0+c*x_0+c*f*x_0*x_0-d*x_0*tf;%tf是脉冲项，这个系统实际上是既有时滞，又有脉冲项，因为这
                                                                              %程序已经有脉冲项了，所以就想把它再改成含时滞的程序，就是改不来。
    y_e1 = 1-2*a*x_0-b*y_0*y_0;
   
   
    %e2
     x_h = x_0 + 0.5*h*x_e1;
     y_h = y_0 + 0.5*h*y_e1;
   
     x_e2 =2*b*x_h*y_h+c*x_h+c*f*x_h*x_h-d*x_h*tf;
     y_e2 = 1-2*a*x_h-b*y_h*y_h;
     
     
        %e3
     x_h = x_0 + 0.5*h*x_e2;
     y_h = y_0 + 0.5*h*y_e2;
   
     x_e3 =2*b*x_h*y_h+c*x_h+c*f*x_h*x_h-d*x_h*tf;
     y_e3 = 1-2*a*x_h-b*y_h*y_h;
   
     
        %e4
     x_h = x_0 + h*x_e3;
     y_h = y_0 + h*y_e3;
   
     x_e4 =2*b*x_h*y_h+c*x_h+c*f*x_h*x_h-d*x_h*tf;
     y_e4 =1-2*a*x_h-b*y_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);
   
    x_0 = x_1;
    y_0 = y_1;
   
   
end

X = X(1:end);
Y = Y(1:end);

figure(1);
plot(X(1:end))
grid off
xlabel('t');
ylabel('x_1');

figure(3);
plot(Y(1:end))
grid off
xlabel('t');
ylabel('y_1');