%this routine solves the system of odes for a chemostat with a food chain
function foodchain
clear
global  K1 K2 c1 c2 

%specify parameter values
r = 1;
K1 = 0.3;
c1 = 2.5;
K2 = 0.3;
c2 = 3;


%specify the temporal output points
tspan = [0:.1:100];

%specify initial data

y0 = [0.8,.1,.1];

[T,Y] = ode23(@foodchain_deRHS,tspan, y0);  
figure(1)
plot(Y(:,2),Y(:,3))
axis([0 1 0 2])
figure(2)
plot(T,Y)


function F_prime=foodchain_deRHS(t,y)
global K1 K2 c1 c2 

s = y(1);
u = y(2);
v = y(3);

Sp = (1-s) - c1*s*u/(K1+s);
Up = c1*s*u/(K1+s) - c2*u*v/(K2+u)-u;
Vp =  c2*u*v/(K2+u)-v;

F_prime = [Sp;Up;Vp];
