% This is the main function to solve 
% dN/dt = k C N in p.119, Mathematical Models in Biology
% dC/dt = -alpha k C N

global k alpha

k = 1; 
alpha = 0.5;

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

t0 = 0;
tfinal = 3;
N_ini = 5;
C_ini = 4;
y_ini = [N_ini C_ini];

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



% plot N at different T
plot(T, y(:,1),T,y(:,2));
legend('N','C')
xlabel t


% plot exact solution
C0 = C_ini+alpha*N_ini;
r = k*C0;
B = C0/alpha;
Nexact = N_ini*B./(N_ini+(B-N_ini)*exp(-r*T));
Cexact = -alpha*Nexact+C0;

figure(2)
plot(T,Nexact,'r',T,Cexact,'r')
legend('N exact solution','C exact solution')
xlabel t, ylabel('exact solution of N')