求助：如何解带时滞的微分方程组

kevin2008

请教高手指点，如何解如下带时滞的微分方程组：

ds(t)/dt=δ*N2(t)-β(t)*s(t)*i(t)+α*i(t)-μ*s(t-T)*j(t-T)-d1*s(t)
di(t)/dt=(t)*s(t)*i(t)-λ*i(t-T)-α*i(t)-d2*i(t)
dr(t-T)/dt=λ*i(t-T)-d3*r(t-T)
dq(t-T)/dt=μ*s(t-T)*j(t-T)-d4*q(t-T)

其中，β(t)=β0*[1-i(t)/N1(t)]^η
     j(t)=i(t)+r(t)
     N1(t)=q(t)+s(t)+i(t)+r(t)
     N=N1(t)+N2(t)

在仿真中取参数N1(0)==300000， =0.1/N1(0)，μ=0.04/N1(0)，η=2，i(0)=1，s(0)=N1(0)-i(0)，λ=0.006，δ=0.003，α=0.001，d1=d2=d3=d4=0.0005，t<T时，r(t)=q(t)=0;

我是这样解得：
ffdu.m
function dy=ffdu(t,y,Z)
N=6*10^5;
N1=y(1)+y(2)+y(3)+y(4);
N2=N-N1;
beta0=0.1/(3*10^5);
mu=0.04/(3*10^5);
lambda=0.006;
delta=0.003;
alpha=0.001;
d1=d2=d3=d4=0.0005;
eta=2;
beta=beta0*(1-y(2)/N1)^eta;

dy=[delta*N2-beta*y(1)*y(2)+alpha*y(2)-mu*Z(1)*(Z(2)+Z(3))-d1*y(1);...
    beta*y(1)*y(2)-lambda*Z(2)-alpha*y(2)-d2*y(2);...
    lambda*Z(2)-d3*Z(3);...
    mu*Z(1)*(Z(2)+Z(3))-d4*Z(4)];

main.m
clear
tt=0:1:800;
options=ddeset('RelTol',1e-4,'AbsTol',1e-7);
sol=dde23('ffdu',[50;50;50;50],[3*10^5-1;1;0;0],tt,options);
plot(sol.x,sol.y(1,:),':k',sol.x,sol.y(2,:),'-r',sol.x,sol.y(3,:),'-.m');
axis([0 800 0 4*10^5]);

提示错误：
??? Error: File: D:\matlab\work\ffdu.m Line: 10 Column: 6
The expression to the left of the equals sign is not a valid target for an assignment.

Error in ==> dde23 at 224
f0 = feval(ddefun,t0,y0,Z0,varargin{:});

Error in ==> du at 5
sol=dde23('ffdu',[50;50;50;50],[3*10^5-1;1;0;0],tt,options);

请高手指点一下，我的解法存在什么问题，多谢！