Fundamental solution of Heat equation For every bounded function $u_0\in C(\mathbb R)$ the function $$u(t,x)=\frac{1}{\sqrt{4\pi t}}\int_{-\infty}^{\infty}e^{-\frac{-(x-y)^2}{4t}}u_0(y)dy$$ is a solution of $\dfrac{\partial v}{\partial t}(t,x)=\dfrac{\partial^2v}{\partial x^2}(t,x)$
and $\lim\limits_{t \to 0}u(t,x)=u_0(x)$.
To prove that is it just necessary to differentiate $$\frac{1}{\sqrt{4\pi t}}e^{-\frac{-(x-y)^2}{4t}}u_0(y)\;?$$ Because when I differentiate it one time to $t$ and two times to $x$ the results are not the same (they differ in the sign). Is there a better method to prove it is a solution?
That $\lim\limits_{t \to 0}u(t,x)=u_0(x)$ is from the physical point of view clear, from the mathematical not when I let $t \to 0$ in the integral term. 
 A: Differentiating the expression for $u(t,x)$ should work, just write it out again more carefully and you should find your mysterious negative sign.  You can also do this "symbolically" by letting $H(t,x)=\frac{1}{\sqrt{4\pi t}}e^{-x^2/4t}$, verifying that $H$ solves $u_t=u_{xx}$ (by taking derivatives carefully) and writing the integral for $u(t,x)$ as: 
$$
u(t,x)=\int_{-\infty}^\infty H(t,x-y)u_0(y)dy
$$
so that 
$$
u_t=\int H_t(t,x-y)u_0(y)dy\\
u_{xx}=\int H_{xx}(t,x-y)u_0(y)dy\\
\Rightarrow u_t-u_{xx}=\int (H_t(t,x-y)-H_{xx}(t,x-y))u_0(y)dy=0
$$ since $H_t=H_{xx}$.
As for showing rigorously that $\lim_{t\rightarrow 0}u(t,x)=u_0(x)$, this requires quite a bit more effort. I would consult a book on basic Fourier analysis (like Stein and Shakarchi) or PDE - the solution is too long to write out here.  You have to make a change of variables, then break the integral up into different parts and carefully control the size of each part in order to show that it goes to $u_0$.  Essentially what you're doing is showing that this function $H$ (the "heat kernel") is an approximation to the identity, or in other words, $\lim_{t\rightarrow 0}H(t,x)=\delta_0(x)$, the Dirac delta distribution. 
A: It was requested by another user that I post the method to show that $$\lim_{t \to 0}u(t,x)=u_0(x)$$ assuming that $u_0$ has a derivative, not necessarily continuous, which is bounded on every finite interval.  I'll also asume that
$$
\int_{-\infty}^{\infty} |u_0(x)|\,dx < \infty.
$$
My apologies if you don't find this relevant.
Fix $x$.  Let $\lambda = 1/(4t)$ and make the change of variables $s = y-x$.  For convenience, define $v_x(s) = u_0(x+s)$.  This gives
$$
\begin{align*}
\int_{-\infty}^{\infty} e^{-(x-y)^2/(4t)} u_0(y)\,dy &= \int_{-\infty}^{\infty} e^{-\lambda s^2} v_x(s) \,ds \\
&= \int_{-\infty}^{\infty} e^{-\lambda s^2} \Bigl[v_x^{\text{odd}}(s) + v_x^{\text{even}}(s)\Bigr] \,ds \\
&= 2\int_{0}^{\infty} e^{-\lambda s^2} v_x^{\text{even}}(s) \,ds,
\tag{1}
\end{align*}
$$
where
$$
v_x^{\text{even}}(s) \triangleq \frac{v_x(s)+v_x(-s)}{2} \qquad \text{and} \qquad v_x^{\text{odd}}(s) \triangleq \frac{v_x(s)-v_x(-s)}{2}.
$$
Making the change of variables $s = \sqrt{r}$ in $(1)$ gives
$$
\begin{align*}
2\int_{0}^{\infty} e^{-\lambda s^2} v_x^{\text{even}}(s) \,ds &= \int_{0}^{\infty} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr \\
&= \int_{0}^{\delta} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr + \int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr,
\qquad \tag{2}
\end{align*}
$$
where $0 < \delta < \infty$ is fixed.  We see that the right integral in $(2)$ is exponentially small as $\lambda \to \infty$ by the calculation
$$
\left|\int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr\right| \leq e^{-\delta \lambda} \int_{\delta}^{\infty} r^{-1/2} \left|v_x^{\text{even}}\left(\sqrt{r}\right)\right| \,dr.
\tag{3}
$$
The finiteness of this last integral follows from the assumption that $u_0$ is integrable.
Now, by assumption we can appeal to Taylor's theorem with remainder to write
$$
\begin{align*}
v_x(s) &= v_x(0) + R_x(s) \\
&= u_0(x) + R_x(s),
\end{align*}
$$
where
$$
|R_x(s)| \leq \sup_{|\tau| < \delta} |u_0'(x+\tau)| \cdot |s|
$$
for $|s| \leq \delta$, so that
$$
v_x^{\text{even}}(s) = u_0(x) + \hat{R}_x(s),
\tag{4}
$$
where
$$
|\hat{R}_x(s)| \leq \sup_{|\tau| < \delta} |u_0'(x+\tau)| \cdot |s|
$$
for $|s| \leq \delta$.  Substituting $(3)$ into the left integral in $(2)$ we get
$$
\int_{0}^{\delta} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr = u_0(x) \int_{0}^{\delta} e^{-\lambda r} r^{-1/2} \,dr + \int_{0}^{\delta} e^{-\lambda r} r^{-1/2} \hat{R}_x\left(\sqrt{r}\right) \,dr.
\quad \tag{5}
$$
We can bound the integral on the right,
$$
\left|\int_{0}^{\delta} e^{-\lambda r} r^{-1/2} \hat{R}_x\left(\sqrt{r}\right) \,dr\right| \leq \sup_{|\tau| < \delta} |u_0'(x-\tau)| \int_0^\delta e^{-\lambda r}\,dr = O(1/\lambda),
\tag{6}
$$
and the integral on the left can be rewritten as
$$
\begin{align*}
\int_{0}^{\delta} e^{-\lambda r} r^{-1/2} \,dr &= \int_{0}^{\infty} e^{-\lambda r} r^{-1/2} \,dr - \int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} \,dr \\
&= \sqrt{\frac{\pi}{\lambda}} - \int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} \,dr,
\tag{7}
\end{align*}
$$
where
$$
0 \leq \int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} \,dr \leq \frac{1}{\sqrt{\delta}\lambda} e^{-\delta\lambda}
\tag{8}
$$
is exponentially small as $\lambda \to \infty$.
Finally, combining the estimates $(3)$, $(6)$, and $(8)$ with the expressions $(1)$, $(2)$, $(5)$, and $(7)$, we get
$$
\int_{-\infty}^{\infty} e^{-\lambda\,(x-y)^2} u_0(y)\,dy = \sqrt{\frac{\pi}{\lambda}}\,u_0(x) + O\left(\frac{1}{\lambda}\right)
$$
as $\lambda \to \infty$.  Equivalently,
$$
\frac{1}{\sqrt{4\pi t}} \int_{-\infty}^{\infty} e^{-(x-y)^2/(4t)} u_0(y)\,dy = u_0(x) + O\left(\sqrt{t}\right)
$$
as $t \to 0^+$.
(The result is surely true for lighter assumptions on $u_0$, perhaps if $|u_0(x)| \leq a e^{b x^2}$.  This could add to the calculations considerably and I wanted to keep things more simple than complicated here.)
