clear all
% set up global variables for use in the DE routine.
global q un uc Kn Kc Vn Vc r
global yv nv uv cv
global icount

t0  = 0;     %initial time
t1  = 40;    % final time
c0  = 0;     % initial carbon conc
n0  = 0;     % initial nitrogen conc
y0  = 4.9;   % initial yeast concentration

% initial state for use in DE solver.

S0  = [c0;n0;y0];

q  = 0.15; % dilution rate
uc = 100;

%
% terSchure data from Figure 1 A.
%
uv = [29 44	61	66	78	90	96	114	118];
yv = [4.9	7	8.1	8.1	7.5	8.1	7.9	8.1	8.2];
nv = [0	0	0	8	21	30	40	55	60];

% terSchure paper observation on carbon.
cv = ones(size(uv))*0.1;

Vn = 53.8607;
Vc = 92.7348;
Kn =  0.1000;
Kc =  4.8231;
r  =  7.4205;

sym = '+o*.xsd^v><ph+o*.xsd^v><ph+o*.xsd^v><ph+o*.xsd^v><ph';

nm = zeros(size(nv));
cm = zeros(size(nv));
ym = zeros(size(nv));

figure(1)

for ii = 1:length(uv)

    un = uv(ii);
    
    [t,St] = ode45('chemostat_dynamics',[t0,t1],S0);
    
    subplot(311)
    plot(t,St(:,1),t,cv(ii)*ones(size(t))),title('carbon')
    
    subplot(312)
    plot(t,St(:,2),t,nv(ii)*ones(size(t))),title('nitrogen')
    
    subplot(313)
    plot(t,St(:,3),t,yv(ii)*ones(size(t))),title('biomass')
    
    pause
    
    cm(ii) = St(end,1);
    nm(ii) = St(end,2);
    ym(ii) = St(end,3);
    
end
figure(2)
plot(uv,nv,uv,nm),title('Nitrogen residual')
figure(3)
plot(uv,cv,uv,cm),title('Carbon residual')
figure(4)
plot(uv,yv,uv,ym),title('Biomass')