Solving a PDE in $1$D I want to solve the below equation
$\partial_t f(x,t)=D \partial_x^2 f(x,t)$
with initial and boundary conditions:
$f(x,0)=\delta(x-x_0)$
$D\partial_x f(x,t)=0$ at $x=L$ and $x=0$;
Can anyone help me with solving it? 
 A: Let $f(x,t) = X(x)T(t)$ then we have (assuming $D\ne 0$)
$$ \frac{T'(t)}{DT(t)} = \frac{X''(x) + \frac{1}{D}X'(x)}{X(x)} = -\lambda $$
where $X(x)$ satisfies the BVP
\begin{cases} 
X'' + \frac{1}{D}X' + \lambda X = 0 \\
X'(0) + \frac{1}{D}X(0) = 0 \\
X'(L) + \frac{1}{D}X(L) = 0
\end{cases}
The characteristic equation is
$$ r^2 + \frac{1}{D}r + \lambda = \left(r+\frac{1}{2D}\right)^2 + \left(\lambda - \frac{1}{4D^2}\right) = 0 $$
Let $\lambda - \frac{1}{4D^2} = \mu^2$, then the solution is
$$ X(x) = e^{-x/2D}\left[A\cos(\mu x) + B\sin(\mu x) \right] $$
Matching up the boundary conditions gives $A=-2D\mu B$ and $\sin(\mu L) = 0$, so we have
$$ X_n(x) = e^{-x/2D}\left[\frac{n\pi}{L} \cos\left(\frac{n\pi}{L}x\right) - \frac{1}{2D}\sin\left(\frac{n\pi}{L}x\right) \right] $$
Along with
$$ T_n(t) = \exp\left(-\left(\frac{n^2\pi^2}{L^2}+\frac{1}{4D^2}\right)Dt\right) $$
Leading to the series solution $u(x,t) = \sum\limits_{n=1}^\infty c_n X_n(x)T_n(t)$. The initial condition requires
$$ u(x,0) = \sum_{n=1}^\infty c_nX_n(x) = \delta(x-x_0) $$
where the constants $c_n$ are given by orthogonality

Edit 3/1: Proof for orthogonality
Here I will show that the eigenfunctions in $X(x)$ are orthogonal with respect to a weighting function. Note that the boundary conditions are equivalent to
$$ X' + \frac{1}{D}X = 0 \implies (e^{x/D}X)' = 0, \quad x = \{0,L\} $$
Let $Y(x) = e^{x/D}X(x)$, then the BVP is reduced to
$$ Y_n'' - \frac{1}{D}Y_n' = -\lambda Y_n, \quad Y_n'(0) = Y_n'(L) = 0$$
or
$$ (e^{-x/D}Y_n')' = -\lambda e^{-x/D}Y_n $$
Following the method similar to this answer we can show that
$$ \int_0^L (e^{-x/D}Y_n')'Y_m dx - \int_0^L (e^{-x/D}Y_m')'Y_n dx = 0 $$
$$ \implies \int_0^L e^{-x/D}Y_nY_m dx = 0 $$
$$ \implies \int_0^L e^{x/D}X_nX_m dx = 0 $$
Now we solve for the constants by
\begin{align} 
\int_0^L e^{x/D}\delta(x-x_0) X_m(x)\ dx &= \sum_{n=1}^\infty c_n \int_0^L e^{x/D} X_n(x)X_m(x)dx \\
\implies e^{x_0/D}X_m(x_0) &= c_m \int_0^L e^{x/D}X_m^2(x) dx \end{align}
Plugging in the expression for $X(x)$ gives
$$ c_m = \frac{\displaystyle e^{-x_0/D}\left[\frac{m\pi}{L} \cos\left(\frac{m\pi}{L}x_0\right) - \frac{1}{2D}\sin\left(\frac{m\pi}{L}x_0\right) \right]}{\displaystyle \int_0^L \left[\frac{m\pi}{L} \cos\left(\frac{m\pi}{L}x\right) - \frac{1}{2D}\sin\left(\frac{m\pi}{L}x\right) \right]^2dx} $$
