clear all
t0 = 0;   %initial time
t1 = 100;   % final time
c0 = 0;  % initial nitrogen conc
N0 = 4.9; % initial biomass
c20 = 0; % initial glucose

S0 = [c0;N0;c20]; % load initial state vector for ode45
%
% the global statement below allows
%the chemostat DE function to know
%all the rate parameters
%
global q u r Kn Kc V u2

% data from the ter Schure paper: feeds, biomass, nitrogen

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];

%
% Model parameters specified
%
q   = 0.15;
r   = 0.14;
V   = 8;
Kn  = 5;
Kc  = 244;
u2  = 100;

ym = zeros(size(yv));

for ii = 1:length(uv)
    u = uv(ii);
    [t,St] = ode45('Chemostat2n',[t0,t1],S0);
    subplot(311)
    plot(t,St(:,1),t,0*t+nv(ii)),title('nitrogen')
    subplot(312)
    plot(t,St(:,3)),title('carbon')
    subplot(313)
    plot(t,St(:,2),t,0*t+yv(ii)),title('biomass')
    pause
    n = length(St(:,2));
    ym(ii) = St(n,2);
end

figure

plot(uv,yv,'ko-',uv,ym,'b')