Pharmacokinetics differential equations; with equal absorption and elimination constants.

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}