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}
Best Answer
Firstly, you can solve the first equation:
$$\frac {dA} {dt} = -k_a \cdot A \Rightarrow A(t) = A_0e^{-k_at},$$
having posed that $A(0) = k_0$. Then, solve the second equation imposing $k_0 = k_a$:
$$\frac {dC} {dt} = - k_a \cdot C + k_a A_0e^{-k_at}.$$
The general solution is $C_g(t) = Me^{-k_at}$, for some $M$. The particular solution should normally be in the form $C_p(t) = Ne^{-k_at}$, too. But since the general solution is similar, then
$$C_p(t) = Nte^{-k_at},$$
and
$$C(t) = Me^{-k_at} + Nte^{-k_at}.$$
Then:
$$C(0) = M = 0,$$
and $$\frac {dC} {dt} = Ne^{-k_at}-Ntk_ae^{-k_at} = -Ntk_ae^{-k_at} + k_aA_0e^{-k_at}.$$
Hence, $N=k_a A_0$ and your unique solution is $C(t) = k_a A_0 te^{-k_at}$
Another way...
Consider the following:
$$\lim_{k_a \to k_0}\frac{A_0 \cdot k_a \cdot (e^{-k_a \cdot t}-e^{-k_0 \cdot t})}{k_a-k_0} = \lim_{k_a \to k_0}\frac{A_0 \cdot k_a \cdot e^{-k_a \cdot t} (1-e^{(k_a-k_0) \cdot t})}{k_a-k_0} =\\= A_0 k_0 t e^{-k_a \cdot t} ,$$
since
$$\lim_{x \to 0} \frac{1 - e^{-xt}}{x} = t$$