function dx = chemostatfunction(t,x)
global q u1 u2 V1 V2 R K1 K2

dx = zeros(size(x));

dx(1) = (q*u1) - (q*x(1)) - (x(3)*V1)*(x(1)/(x(1)+K1))*(x(2)/(x(2)+K2));
dx(2) = (q*u2) - (q*x(1)) - (x(3)*V2)*(x(1)/(x(1)+K1))*(x(2)/(x(2)+K2));
dx(3) = (x(3)*R)*(x(1)/(x(1)+K1))*(x(2)/(x(2)+K2));

end