function dx = chemostatdynamics(t,x)
global q u V K r 

dx = zeros(size(x));

dx(1) = ((q*u)-(q*x(1))-x(2)* (V*x(1))/(x(1)+K));
dx(2) = -(q*x(2))+(r*x(2))*((V*x(1))/(x(1)+K));
dx(3) = ((q*u)-(q*x(3))-x(2)* (V*(x(1)+x(3)))/(x(1)+x(3)+K));

% c=x(1)
% y=x(2)  

end






global  q u V K r                                              
q=3;
u=2;
V=10;
K=30;
r=4; 

                                                                                                                                                                                                                                                            
tt = 0:0.1:10;
x0 = [5;6;1];

[t,x] = ode45('chemostatdynamics',tt,x0);

plot(t,x)
xlabel('time')
ylabel('abundance')
title('Population and food source over time')
legend('Population','Food','Food 2')