龙格库塔法求解延时微分方程matlab

clc;
L=[-10 1 2;3 -12 0;1 2 -9];
M=[-1 0 3;1/2 1/2 1;1/2 3/2 0];
N=1/72*[-1 5 2;4 0 3;-2 4 1];         %生成三个已知矩阵
B=110;                                %计算步数
h=0.1;                                %计算步长
tau=1;                                %时间延迟
m=tau/h;                              %延时时间计算步数
for i=1:m                             %初始化参数以及u
    k1（:i）=[0 0 0];
    k2（:i）=[0 0 0];
    k3（:i）=[0 0 0];
    k4（:i）=[0 0 0];
    u（:i） =[1 1 1];
end
u（:11）=[1 1 1];
for n=11:B                            %计算k1k2k3k4和u在第n步的取值
    k1（:n）=h*L*u（:n）+h*M*u（:n