I am trying to solve the following coupled-diffusion equation \begin{align} \frac{\partial a(x,t)}{\partial t}&= D \frac{\partial^2 a(x,t)}{\partial x^2} \\ \frac{\partial b(x,t)}{\partial t}&=-v \frac{\partial b(x,t)}{\partial x}-\mu a(x,t) b(x,t) \\ \frac{\partial c(x,t)}{\partial t}&= \gamma \frac{\partial^2 c(x,t)}{\partial x^2}+k b(x,t) \end{align}
for the functions $a(x,t), b(x,t)$ and $c(x,t)$, subject to boundary conditions $a(x,t\rightarrow \infty)=0$, $\frac{\partial c(x,t)}{\partial(x)}|_{x=0}=-\frac{\partial a(x,t)}{\partial x}|_{x=0}$,$a(L,t)=a_L$, $b(L,t)=b_L$, $c(L,t)=c_L$, and $\frac{\partial c(x,t)}{\partial x}|_{x=L}=0$, where $D, v, \mu,\gamma, k$, $a_L$, $b_L$ and $c_L$ are constants.
I tried to attack this problem using, first, the separation of variables method. Although I do obtain a "simple" answer for $a(x,t)$ written in terms of exponentials, I am stuck at the second equation. More precisely, this term with $-\mu a(x,t) b(x,t)$ became a nightmare to me. I would try Fourier transform, which could work because I have the ansatz for a$(x,t)$, but I really don't know how to proceed from that.
Here's one approach. I've indicated the methods used, so you can search for them online. A step by step solution would be far too lengthy for this post.
The solution for $A(x,t)$ is
$$ A(x,t)=\int\limits_{-\infty}^{\infty}dx' \ G_0(x,x',t,0)A_0(x') $$
Where $G_0$ is the free space Green's function of the diffusion equation
$$ G_0(x,x',t,t')=\frac{H(t-t')}{\sqrt{4 \pi \gamma (t-t')}} \exp \left[ -\frac{(x-x')^2}{4 \gamma (t-t')}\right] $$
Where $H$ is the Heaviside step function. So that we have
$$ A(x,t)=\frac{1}{\sqrt{4 \pi D t}}\int\limits_{-\infty}^{\infty}dx' \ A_0(x') \exp \left[ -\frac{(x-x')^2}{4Dt}\right] $$
Similar to what was posted in the comments, but note that this is an integral over the initial condition $A_0(x)$. The same solution can be found by Fourier analysis or separation of variables. I'm writing it like this because $G_0$ will be useful later. Fourier transforms on the $B$ equation won't be fruitful because the transform of $AB$ is a convolution. Using the method of characteristics, we write
$$ \partial_tB+\nu\partial_xB=-\mu AB $$
While the chain rule says
$$ \frac{dB}{ds}=\frac{dx}{ds} \partial_x B +\frac{dt}{ds}\partial_t B $$
So the characteristic curves are given by $\frac{dx}{ds}=\nu$, $\frac{dt}{ds}=1$, and $\frac{dB}{ds}=-\mu AB$. With the boundary data on $B(L,t)=B_L$, we parameterize the curves such that $x(s=0)=L$. The first two equations yield
$$ x(s)=\nu s +L \\ t(s)=s+t_0 $$
The equation for $B$ may be integrated
$$ \frac{dB}{B}=-\mu A(s) ds $$
Where $A(s)=A(x(s),t(s))=A(\nu s+L,s+t_0)$, we have
$$ B(s)=\beta(t_0) \exp \left[-\mu \int\limits_0^s ds' A(s') \right] $$
The unknown function $\beta$ is found by applying $B(L,t)=B_L$
$$ B(x,t)=B_L \exp \left[-\mu \int\limits_0^s ds' A(s') \right]_{t_0=t-\frac{x-L}{\nu} \ ; \ \ s=\frac{x-L}{\nu}} $$
The third equation is a diffusion equation with source term $B(x,t)$. I'm considering the problem for $C$ over the finite domain $0<x<L$. The solution is given by an integral over the diffusion equation Green's function $G(x,t)$. Schematically
$$ C(x,t)=\int dx' \int dt' \ G(x,x',t,t') \ B(x',t') $$
However, we must construct $G$ to fulfill the boundary conditions $C(0,t)=0$ and $\partial_x C(L,t)=0$; there are many ways to do this. I think it's simplest using the eigenfunction expansion
$$ G(x,x',t,t')=\sum_{n=1}^\infty \phi_n(x) \phi_n(x') \exp \left[-\gamma(t-t') \left( n \pi + \pi/2\right)^2 \right] \\ \phi_n(x) =\sqrt{\frac{2}{L}} \sin \left( \frac{\pi x}{L} (n+1/2)\right) $$
It's not 'obvious'. Note that $\phi_n$ automatically satisfy $C(0,t)=0$ and $\partial_x C(L,t)=0$. So we can now forget about those conditions. Finally we have
$$ C(x,t)=\int\limits_{0}^L dx' \int\limits_0^t dt' \ k B(x',t') \ G(x,x',t,t') + \int\limits_0^L dx' \ C_0(x') \ G(x,x',t,0) $$
Where $C_0$ is the initial condition on $C$. I think that's as far as we can go without knowing initial conditions. Finally, there is the condition $\partial_x A(0,t)=\partial_x C(0,t)$: if interpreted as a condition on $C$ for known $A$, makes the problem ill posed, as we have already specified the value of $C(0,t)=0$. If interpreted as a condition on $A$, I think we will get a horrible integro-differential equation that $A_0$ must satisfy; just by taking derivatives and equating them.