(matlab算法及编程咨询：http://shop108557885.taobao.com)

    matlab中求时滞常微分方程组有两个函数，一个是dde23函数，求解时滞函数为常数的情形；另一个是ddesd函数，求解时滞函数为一般函数的情形。下边主要讲dde23函数，ddesd可以参考help文档。

    以控制系统中的时滞常微分方程组为例：

http://s4/mw690/8f5303a74d9ad230cee23&690

代码如下：

function y=controsys()
miu=0.1;beta=0.2;omiga=1;xita=0.3;tao=5;%参数的赋值
sc=[0.6,0.2,0.2];%设定初始值
for i=5:5:100
    sol=dde23(@eq,[1 1],@(t) eqhist(t,sc),[i-5 i]);
    %eq是求解的方程组;eqhist是定义的时滞函数的值[i-tao i];
    %[1 1]是S,I的时滞参数，[i-5 i]是求解区间
    sc=sol.y(:,length(sol.y));
    sc(3,:)=sc(1,:)*xita+sc(3,:);sc(1,:)=sc(1,:)*(1-xita);sc=sc';
    %求解新的初始值
    plot(sol.x,sol.y')
    hold on%画出函数图象
    plot([i i],[sol.y(:,length(sol.y)) sc'])
    hold on%连接跳跃点
end%求解跳跃点
%%
function s=eqhist(t,sc)
s=sc;
%%
function dy=eq(t,y,Z)
miu=0.1;beta=0.2;omiga=1;xita=0.3;tao=5;
dy=[0.6 0.2 0.2];
Slag=Z(1);
Ilag=Z(2);
dy=[miu-beta*y(1).*y(2)-miu*y(1);
    beta*y(1)*y(2)-beta*exp(-miu*omiga)*Slag*Ilag-miu*y(2);
    beta*exp(-miu*omiga)*Slag*Ilag-miu*y(3)];