global q u1 u2 V1 V2 R K1 K2 q = 0.15; u1 = 29; u2 = 100; V1 = 110; V2 = 180; K1 = 2.5; K2 = 10; R = 15; tt = 0:0.1:10; x0= [1,0,0]; [t,x] = ode45('chemostatfunction',tt,x0);