Pharmacokinetics differential equations; with equal absorption and elimination constants. There are three closely related questions for the single compartment model applied to oral dosing in pharmacokinetics:


*

*The time evolution of blood concentrations, $t=0$ to $t=\infty$, given one and only one dose to the stomach at $t=0$; and,

*The equilibrium state time evolution of blood concentrations in between doses, $t=0$ at the instant a dose is taken and $t=\tau$ at the instant the next dose is taken (which then sets $t=0$ again), assuming regularly spaced doses, given that enough time has passed to reach that equilibrium state ($t \to \infty$); and,

*The fuller time evolution of blood concentrations, $t=0$ to $t=\infty$, also assuming regularly spaced doses, and also given that the first dose starts at $t=0$ (each dose occurs at $t=n\tau$, where integer $n \ge 0$ and where $\tau$ is the dose interval.)


In each of the above cases, there is an assumed single absorption rate ($k_a$) to model the transfer from the stomach (mouth/tongue, esophagus, stomach, and the rest of the gut) to the blood stream and another assumed single elimination rate ($k_e$) to model the removal from the blood stream. Let's call $A$ the amount of the dose remaining in the gut and awaiting absorption into the blood and $E$ the amount residing in the blood stream and awaiting elimination. (We'll avoid confounding issues about blood volume and concentration, effectiveness, etc.)

The set up for case 1 above is simple:
\begin{equation}
dA = -k_a \cdot A \, dt \\
dE = k_a\cdot A \, dt - k_e \cdot E \, dt
\end{equation}
and I'm able to solve this when $k_a \ne k_e$ and when $k_a = k_e$.
Just to dot the i, regarding case 1, here's the logic I've applied. (I apologize for doing the obvious here.) Assume that D is the single dose amount, which occurs at $t=0$. So the initial condition is $A_0 = D$.
For the stomach:
\begin{equation}
dA = -k_a \cdot A \, dt, \,\,\,\text{where $A_0=D$} \\
\frac{dA}{A} = -k_a \, dt \\
\int \frac{dA}{A} = \int -k_a \, dt \\
ln(A_t) = -k_a \cdot t + C_0 \\
A_t = D \cdot e^{-k_a \cdot t}
\end{equation}
For the blood:
\begin{equation}
dE = k_a\cdot A \, dt - k_e \cdot E \, dt, \,\,\,\text{where $E_0=0$} \\
dE = k_a\cdot D \cdot e^{-k_a \cdot t} \, dt - k_e \cdot E \, dt \\
\frac{dE}{dt} = k_a\cdot D \cdot e^{-k_a \cdot t} - k_e \cdot E \\
\frac{dE}{dt} + k_e \cdot E = k_a\cdot D \cdot e^{-k_a \cdot t} \\
\text{setting integrating factor $\mu=e^{\int k_e \, dt}=e^{k_e \cdot t}$}, then \\
E = \frac{1}{\mu} \int_0^t \mu \cdot k_a \cdot D \cdot e^{-k_a \cdot s} \, ds \\
E = D\cdot e^{-k_e \cdot t}\cdot k_a \int_0^t e^{\left(k_e-k_a\right) \cdot s} \, ds \\
\text{which resolves one of two ways,} \\
E_t = D \cdot k \cdot e^{-k_e \cdot t} \cdot t, \,\,\,\text{where $k = k_a = k_e$} \\
E_t = D \cdot \frac{k_a}{k_a - k_e} \cdot \left( e^{-k_e \cdot t} - e^{-k_a \cdot t} \right), \,\,\,\text{where $k_a \ne k_e$}
\end{equation}
I get that far without difficulty.
Per the above development, I've also found two related questions have already been asked and answered on stackexchange. These are the relatively clearly stated single compartment model, single dose, absorption rate equals elimination rate and the somewhat more confused question on the same single compartment model, single dose, absorption and elimination rates not necessarily the same. Unfortunately, neither of these address case 2 and case 3. Since I already understand case 1, those cases on stackexchange don't further my understanding at all.

I'm particularly interested in developing case 3 above for the situation where $k_a = k_e$. (I have the book answer, but not solution method, when $k_a \ne k_e$, so if I had those solution steps I could probably intercede and find the other case.) I tap out approximately where I'm able to find the equilibrium equation over time for the stomach value, $A$, for case 2, but not the equilibrium equation for $E$ (though I know what it should be, again from stock answers I've already found.)
I also understand there isn't a single approach, but actually several areas in mathematics which can be applied depending upon where my comfort level is at. For these purposes, assume ONE year of undergrad mathematics: basically MTH 251, 252, and 253. This means NO serious use of Laplace transforms (though a glancing familiarity with the idea), no familiarity with variations of parameters methods, and also very little experience with the application of matrix methods here (though I understand them as applied to solving combinations of finite linear equations.)
I need help both with the setup and the solution process for case 3 but case 2 might be used to prepare for it or else as a check to see if it falls out of the solution for case 3 where $t \to \infty$. (Case 1, of course, should match up with case 3's equation where $0 \le t < \tau$.)
I apologize in advance for any perceived lack of effort or clarity in posing this question. Please accept my assurance that I've spent dozens of hours testing my own skills and attempting to find an appropriate resource that I could learn from before posting here. I'm now at the point of hoping there is someone interested enough to help educate me about solving this case-3 problem.

For those needing/wanting to dig into my motivation:
This pharmacokinetics question actually develops because my daughter suffers from grand mal seizures and takes drugs to help manage them. This is my personal interest here and has nothing whatever to do with homework! (It's been several decades since my first year calculus coursework.) I've taken the time to read what I could over many days' worth of searches (and hand to paper work as well, of course) and still find that I'm out of my depth in being able to develop the solutions myself.
That's frustrating because, while I can find some answers stated in tables for certain questions, I can't create my own solutions to related questions as I simply lack the insight either in properly forming the problems or else find I'm lacking the methods in solving them.

It may help that I have found elsewhere (without knowing how it is achieved) that a solution for case 3, where $k_a \ne k_e$, is:
\begin{equation}
E_t = D \cdot \frac{k_a}{k_a - k_e} \cdot \left[ \left( \frac{1 - e^{-n \cdot k_e \cdot \tau}}{1 - e^{-k_e \cdot \tau}} \right) \cdot e^{-k_e \cdot t} - \left( \frac{1 - e^{-n \cdot k_a \cdot \tau}}{1 - e^{-k_a \cdot \tau}} \right) \cdot e^{-k_a \cdot t} \right]
\end{equation}
where $n$ is the number of doses.
In the case where $n=1$, you can see how this nicely resolves into:
\begin{equation}
E_t = D \cdot \frac{k_a}{k_a - k_e} \cdot \left( e^{-k_e \cdot t} - e^{-k_a \cdot t} \right)
\end{equation}
which is the solution for case 1 where $k_a \ne k_e$ (but not where $k_a = k_e$, where it is instead $E_t = D \cdot t \cdot e^{-k \cdot t}$, where $k = k_a = k_e$)
In the case where $n \to \infty$, you can see how this nicely resolves into:
\begin{equation}
E_t = D \cdot \frac{k_a}{k_a - k_e} \cdot \left[ \frac{e^{-k_e \cdot t}}{1 - e^{-k_e \cdot \tau}}  - \frac{e^{-k_a \cdot t}}{1 - e^{-k_a \cdot \tau}}  \right]
\end{equation}
which is the solution for case 2 where $k_a \ne k_e$ (but not where $k_a = k_e$.)
Note that this case 3 situation (where I can't develop my own work, sadly) only works where $k_a \ne k_e$. I'd like to find the answer where $k_a = k_e$. But I'd also, of course, like to know how to solve it either way. Not just the answers, but the approach.
 A: The (slightly modified) expression that you already have,
represents the concentration of the drug 
in plasma after the $n$-th dose is given 
(that is, parameter $t$ starts from zero 
after $(n-1)$ intervals $\tau$ passed),
\begin{align} 
C_n(t) &= 
\frac{\zeta_0}{k_a - k_e} 
\cdot \left[ \left( \frac{1 - \exp(-k_e n\tau)}
{1 - \exp(-k_e \tau)} \right) 
\cdot \exp(-k_e t)\right. 
\\
&-\left.C_n(i+1,\tau,k_a,\zeta_0)(
 \left( \frac{1-\exp(-k_a n\tau)}
{1 - \exp(-k_a  \tau)} \right) 
\cdot \exp(-k_a t) 
\right]
\tag{1}\label{1}
.
\end{align}
Here $\zeta_0$ is the universal model constant, 
the initial slope of the single-dose concentration curve.
\begin{align} 
\zeta_0&=C_0\cdot k_a=\frac{D}{V}\cdot k_a
,
\end{align}
It encapsulates the tricky parameters
like the apparent volume of distribution $V$,
the apparent initial concentration $C_0$
as well as provides a symmetry 
to the expression and invariance 
to possible flip-flop condition
(when it happens that $k_a<k_e$).
So, if we define a function $f$ as
\begin{align} 
f(x)&=
\zeta_0
\cdot \left( \frac{1 - \exp(-x n\tau)}
{1 - \exp(-x \tau)} \right) 
\cdot \exp(-x t) 
,
\end{align}
then expression \eqref{1} in terms of $f$ is
\begin{align} 
\frac{f(k_e)-f(k_a)}{k_a-k_e}
&=
-\frac{f(k_a)-f(k_e)}{k_a-k_e}
.
\end{align}
In order to handle the case when $k_a=k_e$, we just need 
to find a limit
\begin{align} 
C_n(t)|_{k_e=k_a}&=
-\lim_{k_e\to k_a}\frac{f(k_a)-f(k_e)}{k_a-k_e}
\tag{2}\label{2}
.
\end{align}
Cn(i+1,tau,ka,zeta0)(
But \eqref{2} is just a definition of the derivative of $f$, hence
\begin{align}
C_n(t)|_{k_e=k_a}
&=-f'(x)|_{x=k_a}
,
\end{align} 
which we can find to be
\begin{align} 
C_n(t)|_{k_e=k_a}&=
\frac{\zeta_0\exp(-k_a t)}{1-\exp(-k_a\tau)}\cdot
\left[
\frac{(1-\exp(-k_a n\tau))\tau\exp(-k_a\tau)}{1-\exp(-k_a\tau)}
+
t-(n\tau+t)\exp(-k_a n\tau )
\right]
\tag{3}\label{3}
.
\end{align} 
It's easy to check that for $n=1$
expression \eqref{3} gives 
\begin{align}
C_1(t)|_{k_e=k_a}&=
\zeta_0t\exp(-k_at)
,
\end{align}
as expected.
An illustration for $\tau=3$, $k_a=k_e=0.9$, $\zeta_0=1.3$,
$n=1,\dots,6$:

As for the origin of \eqref{1},
it's just a sum of of geometric progression:
\begin{align}
\left( \frac{1 - \exp(-k_e n\tau)}
{1 - \exp(-k_e \tau)} \right) 
\cdot \exp(-k_e t)
&=\sum_{m=0}^{n-1}
\exp(-k_e m\tau)\exp(-k_e t)
\\
&=\sum_{m=0}^{n-1}
\exp(-k_e (m\tau+t))
,
\end{align}
that is, it's a sum of 
all $n$ individual single-dose curves
at a point of $(n-1)\tau+t$.
Edit
An alternative way to get the same expression
is to start from single-dose expression
\begin{align}
C_1(t)|_{k_a=k_e}
&=
\zeta_0\,t\exp(-k_a t)
\end{align}
and just calculate the cumulative effect 
of multiple doses as
\begin{align}
C_n(t)|_{k_a=k_e}
&=
\sum_{m=0}^{n-1}
\zeta_0(m\tau+t)\exp(-k_a(m\tau+t))
.
\end{align}
A: Here is an analysis for the third case.
Let $\mathbb{A} =\left\lbrack \begin{array}{cc}-k_a&0\\ k_a& -k_e \end{array} \right\rbrack$,the state vector 
$$X(t) = \left\lbrack \begin{array}{c} A(t)\\ E(t) \end{array}\right\rbrack,$$
 the input matrix $\mathbb{B} = \left\lbrack \begin{array}{c}1 \\0 \end{array} \right\rbrack$, and the input $$u(t) =D \sum_{i=0}^\infty  \delta(t-i \tau)$$  Then the system (I think) is described by the simple linear time-invariant system of ODEs as follows.
$$
\frac{\mathrm{d}}{\mathrm{d}t} X(t) = \mathbb{A} X(t) + \mathbb{B} u(t)  
$$
For a system this small, it actually is possible to write a closed form solution but it would involves using Eigenvalues and Eigenvectors.  I will skip all this and write the solution, assuming zero initial conditions.  
$$A(t) = D \sum_{i=0}^\infty e^{k_a(t-i\tau)} \Theta(t-i\tau)$$
and 
$$E(t) = D \frac{k_a}{k_a-k_e} \sum_{i=0}^\infty \left( e^{k_a(t-i\tau)}-e^{k_e(t-i\tau)} \right) \Theta(t-i\tau),$$
Where $\Theta$ is the Heaviside function (unit step).  I know you are more interested in knowing how to get the solution, and I can attempt to show you, but it does involve the Laplace Transform.  
Edit
I include the work, so others can check it.
Taking the Laplace Transform of the evolution equation we get
$$(sI-\mathbb{A})X(s) = \mathbb{B} U(s)$$
So we get
\begin{align*}
X(s) &= (sI-\mathbb{A})^{-1} \mathbb{B} U(s)\\
     &= D \left\lbrack\begin{array}{c} (s+k_a)^{-1}\\ k_a(s+k_a)^{-1}(s+k_e)^{-1}\end{array} \right\rbrack \sum_{i=0}^\infty e^{i s \tau} 
\end{align*}
And upon taking the inverse Laplace Transform of this, we get our answer.  Note that we resist the temptation to write $\sum_{i=0}^\infty e^{i s\tau}$ as 
$$\frac{e^s\tau}{e^{s\tau}-1}$$
In the case where $k_a=k_e = k$ we get
$$A(t) = D\sum_{i=0}^\infty e^{-k(t-i\tau))} \Theta(t-i\tau)$$
and 
$$E(t) = D \sum_{i=0}^\infty e^{-k(t-i\tau)}k(t-i\tau) \Theta(t-i\tau)$$
