function [out1, out2, out3]=chemostatmodelus(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);

csection=.5;

%x=X(1).*(X(1)>0);%+(1e-4)*(X(1)<=0);
%y=X(2).*(X(2)>0);
%z=X(3).*(X(3)>0);

x=X(1)+(abs(X(1))+1e-8).*(X(1)<=0);
y=abs(X(2));%X(2)+(abs(X(2))+1e-6).*(X(2)<=0);
z=abs(X(3));%X(3)+(abs(X(3))+1e-6).*(X(3)<=0);
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.*(f1-1-y.*f2)/eps;
    y.*(x.*f2-1-z.*f3);
    del*z.*(y.*f3-1)];
%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/csection-y-z); 
             out1 = [sec1];
             out2 = [1];
             out3 = [-1];
        otherwise
            error(['Unknown flag''' flag ''.'])
    end
end

