% Low frequency relaxation oscillations.
zeta=.1;eps=.01;b1=.48;b2=.2;d1=.42;d2=.4;s1=.25;s2=.25;
tend=1200;
%initial=[y(end,1),y(end,2),y(end,3)]; 
initial=[.9982 .0109 .0846]; %for 'blue'.
options=odeset('e','off','RelTol',1e-8,'AbsTol',[1e-16;1e-16;1e-16]);
%
[t,y]=ode15s('fdchain',[0 tend],initial,options,zeta,eps,b1,b2,d1,d2);
%
% Plot the graphs.
%
figure;
plot3(y(:,1),y(:,3),y(:,2));
% Plot the x-nullcline
hold on;
U= meshgrid(0:.05:1);W=meshgrid(0:.05:1);V=W'; 
X=U;Y=(1-U).*(b1+U)./s1;Z=V; 
mesh(X,Z,Y);
% Plot the z-nullcline on the x-nullcline surface. 
U=0:.05:1; 
yfix=d2*b2/((1-d2)*s1); 
xfix=((1-b1)+sqrt((1-b1)^2-4*(s1*yfix-b1)))/2; 
X=(U-U)+xfix; 
Y=(U-U)+yfix; 
Z=U; 
plot3(X,Z,Y,'k--','LineWidth',1.5);
% plot the y-nullcline on the x-nullcline surface.
hold on;
U=d1*b1/(1-d1):.05:1.03; 
X=U; 
Y=(1-U).*(b1+U)./s1; 
Z=(b2+(1-U).*(b1+U)).*(U./(b1+U)-d1)./s2; 
plot3(X,Z,Y,'k--','LineWidth',1.5); 
axis([0 1 0 1 0 2.5]); 
hidden off;
view(50,20);
text(.7,-.2,-.1,'x');
zlabel('y');
text(1.2,.4,-.1,'z');
% 
