M文件exam1e:
function[value,isterminal,direction] = exam1e(z,x,y)
value =exam1f(z,x,y);
isterminal= zeros(3,1);
direction= -ones(3,1);
end
M文件exam1f:
function v = exam1f(z,x,y)
zlag1 = y(:,1);
v = zeros(3,1);
syms ab B r k1 k2 a1 a2 r1 tau E
a=0.012;
b=0.007;
B=0.6;
r=0.1;
k1=0.02;
k2=0.01;
a1=0.5;
a2=0.6;
r1=0.45;
tau=15;
E=0.006;
v(1)=a*z(2)-r*z(1)-a*zlag1(2)*exp(-r*tau);
v(2)=a*zlag1(2)*exp(-r*tau)+a1*z(3)*z(2)/(1+k1*z(3)+k2*z(2))-B*z(2)^2-E*z(2);
v(3)=r1*z(3)-a2*z(3)*z(2)/(1+k1*z(3)+k2*z(2))-b*z(3);
end
>> options = ddeset('Events','exam1e');
sol =dde23('exam1f',15,[0.2; 0.3; 0.5],[0, 1000],options);
plot(sol.x,sol.y);
xlabel('timet');
ylabel('z(t),x(t),y(t)');