function [out1, out2, out3]=chemostatmodel1(t,X, flag, parameters)
% 
% Use 'ode15s' or 'ode23s' to find solutions to the equation.
%   
%
%eps=.01;mu=1;delta=1;a1=.1;a2=.3;a3=.24;

a1=parameters(1);
a2=parameters(2);
a3=parameters(3);
m1=parameters(4);
m2=parameters(5);
m3=parameters(6);
eps=parameters(7);
del=parameters(8);
%m3=2.5;a3=.2;

%eps=parameters(7);
%del=parameters(8);


%a1=0.08;a2=0.23;a3=.24;m1=10;m2=4;m3=3.5;%m3=2.5;a3=.2;

%eps=1;
%del=1;

x=X(1).*(X(1)>0);
y=X(2).*(X(2)>0);
z=X(3).*(X(3)>0);

%b=1-refx-refz;
%w=tan(pi*(b-y)/2/b);
%y=b-2*b*atan(w)/pi;

f1=m1*(1-x-y-z)./(a1+(1-x-y-z));
f2=m2./(a2+x);
f3=m3./(a3+y);

if nargin <3 | isempty(flag) 
  out1=[(x+1e-6*(x<=0)).*(f1-1-y.*f2)/eps;
    (y+1e-6*(y<=0)).*(x.*f2-1-z.*f3);
    del*(z+1e-6*(z<=0)).*(y.*f3-1)];
else
    switch(flag)
        case 'init'
            out1 = [];
            out2 = [];
            out3 = [];
        case 'events' 
             sec1 = (1-x-y-z);
             sec2 = (1-x/.8-y-z); 
             out1 = [sec1 sec2];
             out2 = [1 1];
             out3 = [-1 1];
        otherwise
            error(['Unknown flag''' flag ''.'])
    end
end

