[y]=Chemos(U0, Um, w, So, X1o, X2o, D, a1, a2, m1, m2, T0, Tm);

function [y]=Chemos(U0, Um, w, So, X1o, X2o, D, a1, a2, m1, m2, T0, Tm);
% function Chemostat_model(U0, Um, w, So, X1o, X2o, D, a1, a2, m1, m2, T0, Tm)
%  modeling of controled chemostat described of system ODE by Michaelis-Menten
% function has 13 positive parameters:
%  U0, Um, w - components of control function U = U0 + Um*sin(w*t)
%  So > 0, X1o > 0, X2o > 0 - inital conditions of system ODE
%  D, a1, a2, m1, m2 - parameters of chemostat
%  T0,Tm - define time span [T0 Tm] of modeling
%if nargin < 1
 % U0=1; Um=0.5; w=1; D=1; a1=0.25; a2=1; m1=2; m2=4.5;T0=0; Tm=10;
 %end
tspan = [T0; Tm];
y0 = [So X1o X2o];
%check(U0,Um,w,D,a1,a2,m1,m2);
%options = odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8 1e-8],'OutputFcn',@odephas3);  
options = odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8 1e-8]);
[t y] = ode45(@cheMehMenten,tspan,y0,options,U0,Um,w,D,a1, a2, m1, m2);

%---------------------------------------------------------------------------
figure;
%plot(t,y(:,1),'-',t,y(:,2),'-',t,y(:,3),'-',t,control_U,'--');
plot(t,y(:,1),'-',t,y(:,2),'-',t,y(:,3),'-');
zoom on;
plotedit on;
solution = [ t , y ];
size_of_solution = size(solution);
dt = (Tm - T0)/size_of_solution(1);
delta_t = num2str(dt);
d_t = ';  dt=';
M=max(y);
d = 'Model:D=';
D_chem = num2str(D);
apar1=';a1='; 
a_1 = num2str(a1);
apar2 = ';a2= '; 
a_2 = num2str(a2);
mpar1 = ';m1= ';
m_1 = num2str(m1);
mpar2 = ';m2= '; 
m_2 = num2str(m2);
parchemostat = strcat(d,D_chem,apar1,a_1,apar2,a_2,mpar1,m_1,mpar2,m_2);
u = 'Control:U =';
Uo = num2str(U0);
plus = '+';
um = num2str(Um);
sinus = '*sin('; w_=num2str(w);closeU = '*t)';
control = strcat(u,Uo,plus,um,sinus,w_,closeU,d_t,delta_t);
Smax=(max(y(:,1)));
X1max=max(y(:,2));
X2max=max(y(:,3));
Ycoord=[Smax X1max X2max ];
title('Grafics:Model of chemostat of Michaelis-Menten S(t)-blue-substrat,X1(t)-green,X2(t)-red - populations');
xlabel( 'On plot - MAX of S(t), X1(t), X2(t) and parameters of control & chemostat    t' ); 
ylabel('S(t) X1(t) X2(t)');
grid;
S = 'Smax=';
Smax=num2str(Smax);
S_leg = strcat(S,Smax);
X1max=num2str(X1max);
X1 = 'X1max=';
X1_leg = strcat(X1,X1max);
X2max=num2str(X2max);
X2 = 'X2max=';
X2_leg = strcat(X2,X2max);
legend(S_leg,X1_leg,X2_leg,control,parchemostat,2);
axis( [ 0, Tm, 0, max(Ycoord)+ max(Ycoord)/3] );




% --------------------------------------------------------------------------
function dydt = cheMehMenten(t,y,U0, Um, w, D, a1, a2, m1, m2)
         U=U0+Um*sin(w*t);
dydt = [
        (U-y(1))*D-(m1*y(1)*y(2))/(a1+y(1))-(m2*y(3)*y(1))/(a2+y(1))
        ((m1*y(1))/(a1+y(1))-D)*y(2)
        ((m2* y(1))/(a2+y(1))-D)* y(3)
       ];
%------------------------------------------------------------------------------------------------------

%author Andrey Egin [e-mail], [e-mail], chief professor Kolosov G. E.
% Moscow State Institute of electronic and mathematics
%created feb 2002