outputting an extra variable from ode system in matlab ode45

Question

I have a ODE system with three variables, $B$- bacteria, $C$-total antibiotic concentration and $C_f$- free drug concentration.

${dB\over dt}=rB-C_fB$.
${dC\over dt}=-\delta C_f$ where
$C_f=1/2*[-(P+K_d-C)+\sqrt{(P+K_d-C)^2+4 K_d C}]$

$r, \delta, K_d, P$ are known constant values. So this system describes bacterial killing through an antibiotic. Here the antibiotic killing should happen based on the free drug concentration ($C_f$. As $C$ decays over time the values of $C_f$ can be computed at each time point. I wrote the code for this in Matlab and can I please know whether

1. based on the way I have written, the correct calculated $C_f$ value gets used to kill the bacteria?
2. Is there a way to out put the value of $C_f$ through the function 'odeModel' rather than separately calculating the value of $C_f$ outside by using the ODE solution obtained for $C$?

y0=[1000,20];  %initial conditions

[t,y]=ode45(@(t,y)odeModel(t,y),0:1:1000,y0);  

function s= odeModel(~,y)
    r1=1;%replication rate of bacteria
    nP=1*10^4;%total concentration of protein binding sites (bound and free)
    Kd=1/(1*10^-4);%disassociation constant of the drug protein complex
    delta=0.08;%decay rate of drug

    s=zeros(2,1);
    s(2)=-(delta/2)*(-(nP+Kd-y(2))+((nP+Kd-y(2))^2+4*Kd*y(2))^0.5);%%the total drug concentration
    Cf=(-(nP+Kd-y(2))+((nP+Kd-y(2))^2+4*Kd*y(2))^0.5)/2;%concentration of the free drug.
    s(1)=r1*y(1)- y(1)*Cf ;%rate of change of bacteria
end

Answer 1

Your Cf looks correct.

y0 = [1000;20]; % initial conditions
[t,y] = ode45(@(t,y) odeModel(t,y),0:1:1000,y0);

To get values of Cf you can change odeModel to give it as an output argument, then once you've solved the ode, you can get the values of Cf by running:

[~,Cf] = odeModel(t,y.')

This is the modified odeModel function:

function [s, Cf] = odeModel(~,y)

r1 = 1; % replication rate of bacteria
nP = 1*10^4; % total concentration of protein binding sites (bound and free)
Kd = 1/(1*10^-4); % disassociation constant of the drug protein complex
delta = 0.08; % decay rate of drug

Cf = (-(nP+Kd-y(2,:))+((nP+Kd-y(2,:)).^2+4*Kd*y(2,:)).^0.5)/2; % concentration of the free drug.

s = zeros(size(y));
s(1,:) = (r1-Cf).*y(1,:); % rate of change of bacteria
s(2,:) = -delta*Cf; % the total drug concentration

end