请高手帮忙看看我的程序错在什么地方？帮忙修改，谢谢了。

xueyongzhou

请高手帮忙看看我的程序错在什么地方？帮忙修改，谢谢了。

events.m

function [value,isterminal,direction]=events(t,y)
global a1 b1 d1 a2 b2 d2
value=(a2*y(2)*y(3)/(1.0+b2*y(2)) - d2*y(3));
isterminal=0;
direction=0;

process.m

clear all;
global a1 b1 d1 a2 b2 d2
a1=5.; a2=0.1; b1=3.; b2=2.; d1=0.25;
% d2=0.0125;%bif par
t0=0.;
tf=20000.;
%y0=[0.66399180 0.20106562 12.053017]
y0=[0.85993 0.07378 10.719]
datout=[0,y0];
trans=800.;
%for d2=0.01434:-0.00002:0.0129
for d2=0.01:-0.00002:0.008
    nsol=size(datout);
    y0=datout(nsol(1),2:4);
    options = odeset('Events',@events,'RelTol',1e-7,'AbsTol',[1e-9 1e-9 1e-9]);
    [T,Y,TE,YE,IE]=ode45(@rm,[t0,tf],y0,options);
    npar=size(TE);
    par=d2*ones(npar(1),1);
    format long
    funcout=[T,Y];
    datoutj=[par,YE];
    datouti=datoutj(trans:npar(1),:);
    datout=cat(1,datout,datouti);
end
%save 'funcout.dat' funcout -ASCII
save 'datout.dat' datout -ASCII
disp('ready')


rm.m

function f = rm(t,y)
global a1 b1 d1 a2 b2 d2
f = zeros(3,1);    % a column vector
f(1)=(y(1)*(1.0-y(1)) - a1*y(1)*y(2)/(1.0+b1*y(1)));
f(2)=(a1*y(1)*y(2)/(1.0+b1*y(1)) - d1*y(2) - a2*y(2)*y(3)/(1.0+b2*y(2)));
f(3)=(a2*y(2)*y(3)/(1.0+b2*y(2)) - d2*y(3));