What is the probability density function over time for a 1-D random walk on a line with boundaries? If a single particle sits on an infinite line and undergoes a 1-D random walk, the probability density of its spatio-temporal evolution is captured by a 1-D gaussian distribution.
\begin{align}
P(x,t)&=\frac{1}{\sqrt{4 \pi D t}}e^{-\frac{(x-x_0)^2}{4Dt}}
\end{align}
However, suppose there are impassable boundaries on the line; on one side, or on both sides. Are there any boundary conditions for which there exists a closed form probability density function for how this particle will behave over time? Any references to such solutions would be extraordinarily helpful.
EDIT. Attempting to generalize Emilio's result below for an arbitrary initial particle position $-L/2 < x_0 < L/2$. 
I had to work it out by example. I found the following "images" were required to account for reflections of an off-center particle at position $x_0$: for the first and second reflections on both sides the new gaussians had to be centered on ($-2L+x_0$, $-L-x_0$, $x_0$, $L-x_0$, $2L+x_0$). From the pattern I think the full solution can be expressed, for all integers $n$, as:
\begin{align}
P(x,t)&=\frac{1}{\sqrt{4 \pi D t}}\sum_{n=-\infty}^\infty e^{-\frac{(x-nL-(-1)^nx_0)^2}{4Dt}}
\end{align}
where the old $x_0$ is now defined as $nL+(-1)^{n}x_0$
 A: This can probably be solved by the method of images, depending on your precise formulation of the problem. The main idea would be to place image particles at the initial time at positions given by treating your impassable boundaries as mirrors; this makes the probability flow at the boundary zero.
To give a more precise formulation, suppose your problem is
$$
\frac{\partial P}{\partial t}=D\frac{\partial^2P}{\partial x^2}\text{ under }\frac{\partial P}{\partial x}(-L/2,t)=0=\frac{\partial P}{\partial x}(L/2,t)\text{ and }P(x,0)=\delta(x),
$$
where I've initially placed the particle in the middle of the barriers for simplicity but this can be altered. The solution is then given, by linearity, by your expression, added up for $x_0=nL$ for all integers $n$:
$$
P(x,t)=\frac{1}{\sqrt{4 \pi D t}}\sum_{n=-\infty}^\infty e^{-\frac{(x-nL)^2}{4Dt}}.
$$
This can be solved exactly in terms of Jacobi theta functions, which makes the calculations and graphing a lot faster, but does not necessarily (at a first go) make this easier to work with:
$$
P(x,t)=\frac{1}{L}\vartheta _3\left(\frac{\pi  x}{L},e^{-\frac{4 D \pi ^2 t}{L^2}}\right).
$$
(For asymetrically placed initial particles, you would have two series of gaussians separated by $2L$, so therefore two theta functions.)
I'm not sure this is very useful by itself, but the method of images is very powerful.
A: Your solution is in fact a particular solution of the 1-Dimensional heat equation, $\frac{\partial  P}{\partial t}  = D \frac{\partial^2  P}{\partial x^2}$,with initial condition $P(x,0) = \delta (x-x_0)$
A traditionnal way to solve this equation is to use Fourier Series
See some solutions in 1-D, like homogeneous equations or inhomogeneous equations, or other examples
