Let us consider the heat equation in one dimension:
$$\dfrac{\partial^2 u}{\partial x^2} = \dfrac{\partial u}{\partial t},$$
defined for $x\in \mathbb{R}$ and $t \geq 0$. Suppose we also impose the condition $u(x,0) = f(x)$.
Then by taking the Fourier transform of the equation we can find out that the solution is
$$u(x,t) = f\ast \mathcal{H}_t(x),$$
where $\mathcal{H}_t(x)$ is the Heat kernel defined by
$$\mathcal{H}_t(x) = \dfrac{1}{\sqrt{4\pi t}}e^{-x^2/4t}.$$
Now, by the very definition of $\mathcal{H}_t$, we have that $\mathcal{H}_t$ is an approximation of the identity as $t\to 0$. This means that by a theorem, $u(x,t)\to f(x)$ as $t\to 0$ uniformly in $x$. In other words, given $\epsilon > 0$ there's $\delta > 0$ such that if $0<t <\delta$ then $|u(x,t)-f(x)|<\epsilon$ for each $x\in \mathbb{R}$.
Now using this I want to show that $u$ is continuous on the closure of the upper half plane $\overline{\mathbb{R}_+^2}$, that is, the set $\mathbb{R}\times [0,\infty)$.
The upper half plane $\mathbb{R}^2_+$ is a subset of this set. Since $u$ is $C^2$ there, it is obviously continuous. Now, when we take the closure I really don't know how to proceed.
I mean, this obviously follows from the fact that $u(x,t)\to f(x)$ uniformly in $x$ as $t\to 0$, but I can't formalize this.
How do I show that $u$ is continuous on the closure of the upper half plane?
First of all, you have to assume that $f$ is continuous and satisfies a growth condition at $\infty$. Then $u(x,t)$ converges uniformly to $f(x)$ as $t\to0$ on compact subsets of $\mathbb{R}$.
To simplify things I will assume that $f$ is bounded, say $|f(x)|\le M$ for all $x\in\mathbb{R}$. Now let $x_0\in \mathbb{R}$. Given $\epsilon>0$ there is a $\delta>0$ such that if $|z-x_0|\le\delta$ then $|f(z)-f(x_0)|\le\epsilon$. Take $x$ such that $|x-x_0|\le\delta/2$. If $|y|\le\delta/2$, then $(x-y)-x_0|\le\delta$.
Then $$\begin{align} |u(x,t)-f(x_0)|&\le\frac{1}{\sqrt{4\,\pi\,t}}\int_{-\infty}^{\infty}e^{-y^2/4t}|f(x-y)-f(x_0)|\,dy\\ &\le\frac{\epsilon}{\sqrt{4\,\pi\,t}}\int_{|y|\le\delta/2}e^{-(x-y)^2/4t}\,dy+\frac{2\,M}{\sqrt{4\,\pi\,t}}\int_{|y|>\delta/2}e^{-(x-y)^2/4t}\,dy\\ &\le\epsilon+\frac{2\,M}{\sqrt{\pi}}\int_{|2\sqrt{t}z-x|>\delta/2}e^{-z^2}\,dx, \end{align}$$ and the last integral converges to $0$ as $t\to0$ uniformly on the interval $|x|\le\delta/2$.