function [ rec ] = BarkaiTimeSeries()
% uses very crude ODE solver to simulate dynamics behaviour of integral
% control model
% alpha is varied randomly with time
% matrix of time series is returned: time | alpha | E | Em | Em'

% initial conditions
E = [100 10 0];
% model parameters
VR = 0.1; VB = 1; Kb = 2; 
dt = 0.5;
rec = zeros(length(0:dt:1000), 5);
c = 1;
residual = 1;
alpha = rand();
% loop through time
for t=0:dt:1000
   NE = E;
   % update concentrations
   NE = NE + dt*[ -VR+VB*E(3)/(Kb+E(3)) VR -VB*E(3)/(Kb+E(3)) ];
   temp = NE(2)+NE(3);
   NE(2) = (1-alpha)*temp;
   NE(3) = alpha*temp;
   residual = NE(3)-E(3);
   E = NE;
   % record concentrations
   rec(c,:) = [c*dt alpha E(1) E(2) E(3)];
   c = c+1;
   if mod(c, 100) == 0
       alpha = rand();
   end
end
stableactivity = E(3);    

end

