\begin{equation}
\frac{dc_1}{d\tau}= \alpha I(1-c_{0}) + c_{1} (-K_{F} - K_{D}-K_{N} s_{0}-K_{P}(1-q_{0}))+ c_{0}(-K_{N} s_{1}+K_{P}q_{1}), \nonumber
\end{equation}
\begin{equation}
\frac{ds_1}{d\tau}= \Lambda_{B} P_{C} (c_{1}(1-s_{0})-c_{0} s_{1})- \lambda_{r} (1-q_{0}) s_{0}, \nonumber
\end{equation}
\begin{equation}
\frac{dq_1}{d\tau}= \frac {P_C}{P_Q} K_{P} ((1-q_{0})c_{1}- c_{0} q_{1}) - \gamma \ q_{0}. \nonumber
\end{equation}
And
\begin{equation}
\frac{dc_2}{d\tau}= - \alpha I c_{1}+ c_{2} (-K_{F} - K_{D}-K_{N} s_{0}-K_{P}(1-q_{0}))-K_N(c_{1}s_{1}+c_{0}s_{2})+K_{P}(q_{1}c_{1}+q_{2}c_{0}), \nonumber
\end{equation}
\begin{equation}
\frac{ds_2}{d\tau}= \Lambda_{B} P_C (c_2(1-s_0)-(c_1 s_1+c_0s_2))- \lambda_r (q_1s_0-s_1(1-q_0), \nonumber
\end{equation}
\begin{equation}
\frac{dq_2}{d\tau}= \frac {P_C}{P_Q} K_P ((1-q_0)c_2-(q_1c_1+ c_0 q_2)) - \gamma \ q_1. \nonumber
\end{equation}
For initial conditions
\begin{equation}
c_0(0)= c(0) = 0.0 \nonumber
\end{equation}
\begin{equation}
s_0(0)= s(0) = 0.02 \nonumber \nonumber
\end{equation}
\begin{equation}
q_0(0)=q(0) = 0.0 \nonumber \nonumber
\end{equation}
and all other terms for $c$, $s$ and $q$ are $0$ at t=$0$ after first terms
\begin{equation}
c_i(0)= 0, \ i>0\nonumber
\end{equation}
\begin{equation}
s_i(0)= 0, \ i>0 \nonumber \nonumber
\end{equation}
\begin{equation}
q_i(0)=0, i>0. \nonumber \nonumber
\end{equation}
The vales of parameters are: $k_f= 6.7*10^{7}$
$ k_d= 6.03*10^8$
$ k_n=2.92*10^9$
$ k_p=4.94*10^9$
$ \alpha =1.14437*10^{-3}$
$I=1200$
$ K_F= k_f * 10^{-9}$
$ K_D= k_d * 10^{-9}$
$ K_N= k_n * 10^{-9}$
$ K_P= k_p * 10^{-9}$
$ P_C= 3 * 10^{11}$
$ P_Q= 2.87 * 10^{10}$
$\lambda_b= 0.0087$
$ \lambda_r =835$
$ \gamma =2.74 $
$ \Lambda_B= \lambda_b *10^{-9}$
I want to find $c_1,s_1,q_1$ and $c_2,s_2,q_2$ also want to plot each $c_1,s_1,q_1$ and $c_2,s_2,q_2$ against $t$ separately in matlab. Any one please help me