[proofplan]
We write the expectation as convolution with the Brownian transition density. Smoothness for positive time follows from differentiating the Gaussian kernel under the integral sign, justified by Gaussian derivative bounds and boundedness of $f$. The [heat equation](/page/Heat%20Equation) follows because the Gaussian kernel itself satisfies $\partial_t p=\frac12\partial_{xx}p$. Continuity at $t=0$ is the approximate-identity property of the Gaussian kernel, proved by splitting the standard normal integral into a compact part and a small tail.
[/proofplan]
[step:Represent the expectation using the Brownian transition density]
For $t>0$, [Brownian motion](/page/Brownian%20Motion) started at $x$ has density
\begin{align*}
p(t,x,y)=\frac{1}{\sqrt{2\pi t}}\exp\left(-\frac{(y-x)^2}{2t}\right)
\end{align*}
with respect to $\mathcal L^1$. Since $f$ is bounded and $p(t,x,\cdot)$ integrates to $1$, the random variable $f(W_t)$ is integrable and
\begin{align*}
u(t,x)=\int_{\mathbb R} f(y)p(t,x,y)\,d\mathcal L^1(y),\qquad t>0.
\end{align*}
At $t=0$, $W_0=x$ $\mathbb P^x$-almost surely, so
\begin{align*}
u(0,x)=\mathbb E^x[f(W_0)]=f(x).
\end{align*}
[/step]
[step:Differentiate the Gaussian convolution for positive time]
Define the heat kernel
\begin{align*}
G:(0,\infty)\times\mathbb R&\to(0,\infty)\\
(t,z)&\mapsto \frac{1}{\sqrt{2\pi t}}\exp\left(-\frac{z^2}{2t}\right).
\end{align*}
Then $p(t,x,y)=G(t,y-x)$, and
\begin{align*}
u(t,x)=\int_{\mathbb R}f(y)G(t,y-x)\,d\mathcal L^1(y).
\end{align*}
Let $K=[a,b]\times[-R,R]\subset(0,\infty)\times\mathbb R$ be compact, with $0<a<b$ and $R>0$. For every pair of integers $r,k\ge0$, differentiating the Gaussian formula shows that
\begin{align*}
\partial_t^r\partial_x^kG(t,y-x)
=P_{r,k}\left(t^{-1},y-x\right)G(t,y-x)
\end{align*}
for a polynomial $P_{r,k}$ depending only on $r$ and $k$. Since $t\in[a,b]$ and $x\in[-R,R]$, there exist constants $C_{r,k,K}>0$, $N_{r,k}\in\mathbb N$, and $c_a>0$ such that
\begin{align*}
\left|\partial_t^r\partial_x^kG(t,y-x)\right|
\le C_{r,k,K}(1+|y|)^{N_{r,k}}\exp(-c_a y^2)
\end{align*}
for all $(t,x)\in K$ and all $y\in\mathbb R$. The function on the right belongs to $L^1(\mathbb R,\mathcal B(\mathbb R),\mathcal L^1)$, and $|f(y)|\le \|f\|_\infty$.
The [dominated convergence theorem](/theorems/4) therefore justifies differentiating under the integral sign any finite number of times on $K$. Since $K$ was arbitrary, $u\in C^\infty((0,\infty)\times\mathbb R)$, and
\begin{align*}
\partial_t u(t,x)&=\int_{\mathbb R}f(y)\partial_tG(t,y-x)\,d\mathcal L^1(y),\\
\partial_{xx}u(t,x)&=\int_{\mathbb R}f(y)\partial_{xx}G(t,y-x)\,d\mathcal L^1(y).
\end{align*}
[/step]
[step:Use the heat-kernel identity to obtain the heat equation]
The Gaussian computation gives
\begin{align*}
\partial_tG(t,z)=\frac12\partial_{zz}G(t,z),\qquad t>0,\ z\in\mathbb R.
\end{align*}
Because $z=y-x$, the second derivative in $x$ equals the second derivative in $z$:
\begin{align*}
\partial_{xx}G(t,y-x)=\partial_{zz}G(t,y-x).
\end{align*}
Using the differentiated integral formulas from the previous step,
\begin{align*}
\partial_tu(t,x)
&=\int_{\mathbb R}f(y)\partial_tG(t,y-x)\,d\mathcal L^1(y)\\
&=\frac12\int_{\mathbb R}f(y)\partial_{zz}G(t,y-x)\,d\mathcal L^1(y)\\
&=\frac12\int_{\mathbb R}f(y)\partial_{xx}G(t,y-x)\,d\mathcal L^1(y)\\
&=\frac12\partial_{xx}u(t,x).
\end{align*}
Thus $u$ solves the [heat equation](/page/Heat%20Equation) on $(0,\infty)\times\mathbb R$.
[/step]
[step:Prove continuity at the initial time]
Let
\begin{align*}
\gamma:\mathbb R&\to(0,\infty)\\
z&\mapsto \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^2}{2}\right)
\end{align*}
be the standard normal density. For $t>0$, the substitution $y=x+\sqrt t\,z$ gives
\begin{align*}
u(t,x)=\int_{\mathbb R} f(x+\sqrt t\,z)\gamma(z)\,d\mathcal L^1(z).
\end{align*}
Fix $x_0\in\mathbb R$ and $\varepsilon>0$. Let $M:=\|f\|_\infty$. By continuity of $f$ at $x_0$, choose $\delta>0$ such that
\begin{align*}
|f(y)-f(x_0)|<\varepsilon
\end{align*}
whenever $|y-x_0|<\delta$. If $|x-x_0|<\delta/2$, then
\begin{align*}
|x+\sqrt t\,z-x_0|<\delta
\end{align*}
whenever $|z|<\delta/(2\sqrt t)$. Hence
\begin{align*}
|u(t,x)-f(x_0)|
&\le \int_{\mathbb R}|f(x+\sqrt t\,z)-f(x_0)|\gamma(z)\,d\mathcal L^1(z)\\
&\le \varepsilon\int_{\mathbb R}\gamma(z)\,d\mathcal L^1(z)
+2M\int_{\{z\in\mathbb R: |z|\ge \delta/(2\sqrt t)\}}\gamma(z)\,d\mathcal L^1(z)\\
&=\varepsilon+2M\int_{\{z\in\mathbb R: |z|\ge \delta/(2\sqrt t)\}}\gamma(z)\,d\mathcal L^1(z).
\end{align*}
The Gaussian tail on the right tends to $0$ as $t\to0^+$. Therefore
\begin{align*}
u(t,x)\to f(x_0)
\end{align*}
as $(t,x)\to(0,x_0)$ with $t\ge0$. Since $x_0$ was arbitrary and $u$ is already smooth for $t>0$, we have $u\in C([0,\infty)\times\mathbb R)$.
[/step]