Reaction diffusion equation solution This has been driving me spare for the last week, and I feel pretty bad for not being able to get a solution because (at least on the face of it) it's a pretty simple equation.
I have the following reaction diffusion equation:
$$\frac{\partial M}{\partial t}=d\frac{\partial^2 M}{\partial x^2}-gM$$
With:
$$\frac{\partial M}{\partial x}(0,t)=-h, \quad M(1,t)=0, \quad M(x,0)=0$$
$$0 \le x \le 1, \quad t>0, \quad d,g,h \ge 0$$
I'm searching for an analytical solution. I've tried separation of variables, I've tried basic transformations, I've looked in books, but I just can't find a solution that satisfies the BCs. I've solved it numerically and (for the values I tried, relating to the larger problem I'm working on) it's a pretty boring curve, so I don't expect there to be anything crazy going on in a solution. 
Please, someone put me out of my misery! I'm fully expecting that I'm missing something obvious or making a ridiculous mistake, but I'd like to see what people answer before I put my working up here.
 A: The boundary condition at $x=0$ is not homogeneous. You cannot apply separation of variables directly. Let $M(x,t)=v(x,t)+\phi(x)$, where $\phi(x)$ is the steady state, that is, the unique solution of
$$
d\,\phi''-g\,\phi=0,\quad\phi'(0)=-h,\quad \phi(1)=0.
$$
Then $v$ satisfies the equation
$$
v_t=d\,v_{xx}-g\,v
$$
with boundary condition $v_x(0,t)=0$, $v(1,t)=0$ and initial value $v(x,0)=-\phi(x)$. Finally, let $v=e^{-gt}w$. Then $w$ satisfies the equation
$$
w_t=d\,w_{xx}
$$
with boundary condition $w_x(0,t)=0$, $w(1,t)=0$ and initial value $w(x,0)=-\phi(x)$. You can now use separation of variables to find $w$ and
$$
M=e^{-gt}\,w+\phi.
$$
As $t\to\infty$, $M$ approaches the steady state solution, as it should be.
A: The problem your describing was solved in Eqs. 13, 14, 15, and Appendix A of https://doi.org/10.1103/PhysRevE.99.042136,
see https://arxiv.org/abs/1902.03963 for a preprint version.
That paper studied the electrolyte Seebeck effect--the tendency of different types of ions to separate when subjected to a temperature difference.
Your $M$ corresponds in that paper to the local ionic charge density $q$ between two flat plates held at different temperatures.
I'm following that paper, where your problem is solved through Laplace inversions.
Original problem:
\begin{align}
\partial_t M &= d \partial_x^2 M-g M\,,\\
\partial_x M(0,t)&=-h\,,\\
M(1,t)&=0\,,\\
M(x,0)&=0\,.
\end{align}
Let's clean up the notation first, $M\to f h$, $t\to t'/d$, and $g\to n^2 d$, and drop primes,
\begin{align}
\partial_t f &= \partial_x^2 f-n^2 f\,,\\
\partial_x f(0,t)&=-1\,,\\
f(1,t)&=0\,,\\
f(x,0)&=0\,.
\end{align}
A Laplace transformation $\hat{f}(x,s)=\mathcal{L}\left\{f(x,t)\right\}=\int_{0}^{\infty}{\rm d}t\, e^{st}f(x,t)$ takes the above equations to
\begin{align}
s\hat{f}(x,s) -f(0,x)&=  \partial_x^2 \hat{f}(x,s)-n^2 \hat{f}(x,s)\nonumber\\
\partial_x^2 \hat{f}(x,s) &=  m^2\hat{f}(x,s)\,,
\end{align}
where we used that $f(0,x)=0$ and where $m^2=n^2+s$.
The boundary conditions turn into
\begin{align}
\partial_x \hat{f}(0,s)&=-\frac{1}{s}\,,\\
\hat{f}(1,s)&=0\,.
\end{align}
The above ODE for $\hat{f}(x,s)$ is solved by
\begin{align}
\hat{f}(x,s)=-\frac{1}{sm}\frac{\sinh[m(x-1)]}{\cosh m}\,.
\end{align}
Determining $f(x,t)=\mathcal{L}^{-1}\left\{\hat{f}(x,s)\right\}$ requires performing an inverse Laplace transformation, which by the residue theorem amounts to
\begin{align}
f(x,t)&=\sum_{s\in s_{\ell}}\text{Res}\left(\hat{f}(x,s)\exp(st),s_{\ell}\right)\,,
\end{align}
where the poles $s_{\ell}=\{s_{0}, s_{g}, s^{\star}_{j} \}$ of $\hat{f}(x, s)$ are located at $s_{0}=0$, $s_{n}=-n^2$, and $s^{\star}_{j}=(m^{\star}_{j})^2 -n^2$ where  $m^{\star}_{j}=\pm i (j-1/2)\pi\equiv\pm i\mathcal{N}_{j}$ with $j=1,2,3\ldots$ (such that $\cosh m^{\star}_{j} =0$).
The pole $s_{0}=0$ gives the steady-state solution,
\begin{align}
\text{Res}\left(f(x,s)\exp(s t),0\right)=-\frac{1}{n}\frac{\sinh[n(x-1)]}{\cosh n}.
\end{align}
For the residue of the pole at $s_{n}=-n^2$, we consider
\begin{align}\label{eq:fkd}
\hat{f}(x,s)\overset{s\to-n^2}=&-\frac{1}{s}\frac{(x-1)}{1}=\frac{1}{n^2}(x-1)\,.
\end{align}
This implies
\begin{align}
\text{Res}\left(\hat{f}\exp(st),s=-n^2\right)=\lim_{s\to-n^2}(s+n^2)\hat{f}(x,s)e^{-n^2t}=0.
\end{align}
For the poles at $s^{\star}_{j}$ we expand
\begin{align}
\cosh (m)\overset{s\to s^{\star}_{j}}=& \frac{\sinh m}{2m}\bigg|_{s=s^{\star}_{j}}\left(s- s^{\star}_{j}\right)\\
\Rightarrow \frac{1}{\cosh (m)}\overset{s\to s^{\star}_{j}}=& \frac{2i(-1)^{j}m^{\star}}{s-s^{\star}_{j}}\,,
\end{align}
where, going to the second line we used $m(s^{\star}_{j})=\pm m^{\star}_{j}$, and $\sinh m^{\star}_{j}=i(-1)^{j+1}$.
We find
\begin{align}
\sum_{j\ge1}\text{Res}\left(\hat{f}\exp(st),s^{\star}_{j}\right)
&=-\sum_{j\ge1}
\text{Res}\left(\frac{2i(-1)^{j}\sinh [m^{\star}(x-1) ] }{ s^{\star}}\frac{\exp(s t)}{s- s^{\star}},s^{\star}_{j}\right)\\
&=-2\sum_{j\ge1}\frac{(-1)^{j}\sin[\mathcal{N}_{j}(x-1)]}{n^2+\mathcal{N}_{j}^{2}}\exp{\left[-\left( n^2+\mathcal{N}_{j}^{2}\right)t\,\right]}.
\end{align}
Collecting the above terms,  we find
\begin{align}
f(x,t)&=-\frac{1}{n}\frac{\sinh[n(x-1)]}{\cosh n}-2\sum_{j\ge1}\frac{(-1)^{j}\sin[\mathcal{N}_{j}(x-1)]}{n^2+\mathcal{N}_{j}^{2}}\exp{\left[-\left( n^2+\mathcal{N}_{j}^{2}\right)t\,\right]},
\end{align}
where, again $\mathcal{N}_{j}=(j-1/2)\pi$.
Let's plot this for $n=10$ and $t=5\cdot10^{-3}$ and $t=1$, and compare to a numerical Laplace inversion of $\hat{f}(x,s)$ for the same times.

I produced this with the following code:
import numpy as np
import matplotlib.pyplot as plt
from mpmath import *

def f(t,x,n):
    fs0 = - np.sinh(n*(x-1))/(n*np.cosh(n))
    fsj=0
    for j in range(1,1000):
        Nj=(j-1/2)*np.pi
        sj=-Nj**2-n**2
        fsj+=  2*(-1)**j*np.sin(Nj*(x-1))*np.exp(sj*t)/sj
    return fs0 + fsj

def F(s,x,n):
    m = mp.sqrt(n**2+s)
    return -mp.sinh(m*(x-1))/(s*m*mp.cosh(m))

mp.dps = 15
def fnum(t,x,n):
    fp = lambda s: F(s,x,n)
    return invertlaplace(fp,t,method='talbot')

x  = np.linspace(0,1,100)
x2 = np.linspace(0,1,20)

n=10
fig, ax = plt.subplots()
ax.plot(x, [f(5e-3,i,n) for i in x],    'r',lw=2,label='$t=5\cdot10^{-3}$, analytical, 1000 terms')
ax.plot(x, [f(1,i,n) for i in x],       'k',lw=2,label='$t=1$')
ax.plot(x2,[fnum(5e-3,i,n) for i in x2],'r',lw=0,marker='o',label='numerical')
ax.plot(x2,[fnum(1,i,n) for i in x2],   'k',lw=0,marker='o',)
ax.set(xlabel=r'$x$', ylabel=r'$f(x,t)$') 
ax.legend()

Returning to the original notation,
\begin{align}
\frac{M(x,t)}{h}&=-\frac{1}{\sqrt{g/d} }\frac{\sinh [\sqrt{g/d} (x-1)]}{\cosh  \sqrt{g/d} }-2\sum_{j\ge1}\frac{(-1)^{j}\sin\!\left[\mathcal{N}_{j}(x-1)\right]}{g/d+\mathcal{N}_{j}^2}\exp{\left[-t\left(g+d\mathcal{N}_{j}^2\right]\,\right]}
\end{align}
