clear all
% set up global variables for use in the DE routine.
global q uy uz Vmax Ky Kz eps delta
global yv nv uv cv
global icount

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

% initial state for use in DE solver.

S0  = [x0;y0;z0];

q   = 0.15; % dilution rate
eps = q/1.1;

%
% 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;
%
% Initial guesses of possible parameter values
%
Vmax  = 10;
Ky  = 0.1;
Kz  = 5;
uz  = 100;
delta = eps/2;

% stack initial guesses into a guess vector

x0 = [Vmax;Ky;Kz;delta;uz];

icount = 0;
options = optimset('MaxFunEvals',1e6,'MaxIter',1e6,'TolFun',1e-4,'TolX',1e-4);
%
% fmincon is the old optimization function.  fminsearch is simpler to use.
%
% x1 = fmincon('simple_fit_fn',x0,[],[],[],[],[1e-4;1e-4;1e-4;1e-4;1e-4],[1e5;1e5;1e5;1e5;1e5],[],options)
x1 = fminsearch('simple_fit_fn',x0);

%
% x1 is the result of the optimization.  Unpack the vector into the
% original variable names
%
Vmax = x1(1);
Ky = x1(2);
Kz = x1(3);
delta = x1(4);
uz = x1(5);

xm = zeros(size(nv));
ym = zeros(size(nv));
zm = zeros(size(nv));

for ii = 1:length(uv)

    uy = uv(ii); % use the feed rate from the list 
    
    [t,St] = ode45('chemostat_dynamics_2n',[t0,t1],S0); % solve the DE
    
    figure % generate a new blank figure window
    
    subplot(311)
    plot(t,St(:,1),t,yv(ii)*ones(size(t))),title('biomass')

    subplot(312)
    plot(t,St(:,2),t,nv(ii)*ones(size(t))),title('nitrogen')

    subplot(313)
    plot(t,St(:,3)),title('??')
    
    pause(0.5)
    
    xm(ii) = St(end,1);
    ym(ii) = St(end,2);
    zm(ii) = St(end,3);
    
end
figure 
plot(uv,nv,uv,ym),title('Nitrogen residual')
figure
plot(uv,cv,uv,zm),title('?? residual')
figure
plot(uv,yv,uv,xm),title('Biomass')