% This is the main function to solve 
% the chem kinetics in p.274



global k1 km1 k2

k1 = 0.1;
km1 = 1.e-5;
k2 = 1;




options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1.e-4 1.e-4]);

t0 = 0;
tfinal =60;


c_ini = [1 2 3 4 5 6];
x0_ini = [1 1 1 1 1 1];
x1_ini = [1 1 1 1 1 1];
p_ini = [1 1 1 1 1 1];


for ind = 1:6

y_ini = [c_ini(ind) x0_ini(ind) x1_ini(ind) p_ini(ind)];

[T,y] = ode45(@chem_kinetic_eqn,[t0 tfinal],y_ini,options);





% plot S, I, Sp, Ip versus T

figure(1);hold on; plot(T,y(:,1),T,y(:,2),T,y(:,3),T,y(:,4));
xlabel('t')
legend('c','x0','x1','p')

figure(2); plot(T, y(:,2)+y(:,3))
ylabel('x_0+x_1')
xlabel('t')
title('conservation of the total occupied and unoccupied receptors')

Kmax = k2*(y(1,2)+y(1,3));
kn = (km1+k2)/k1;

c = y(:,1);
x0 = y(:,2);
x1 = y(:,3);
m_dc_dt = -(-k1*c.*x0+km1*x1);
col = 'bgrcmy'
figure(3); hold on; plot(c, m_dc_dt,col(ind))
xlabel('c')
ylabel('- dc / dt')

pause(3)


end

c = 0:0.1:10;
exact_m_dc_dt =  Kmax*c./(kn+c)
figure(3); hold on; plot(c,exact_m_dc_dt)
xlabel('c')
ylabel('- dc / dt')