How can I prove that the wave function remains normalized as time goes? Exploiting Schrödinger equation and its conjugate we can show that
$$
\dot{\Psi} = \frac{i \hbar}{2m} \nabla^2 \Psi - \frac{i}{\hbar} U \Psi
$$
$$
\dot{\Psi}^* = -\frac{i \hbar}{2m} \nabla^2 \Psi^* + \frac{i}{\hbar} U \Psi^*
$$
So
$$
\frac{\partial |\Psi|^2}{\partial t} = \frac{i \hbar}{2m} (\Psi^* \nabla^2 \Psi - \Psi \nabla^2 \Psi^*)
$$
Exploiting the vector identity $f \nabla^2 g = \nabla \cdot (f \nabla g) - \nabla g \cdot \nabla f$ (with $f$ and $g$ scalar fields) and the divergence theorem (we integrate over $V$), we can rearrange in this way
$$
\frac{d}{dt} \int_V |\Psi|^2 dV = \frac{i \hbar}{2m} \oint_S [\Psi^* \nabla \Psi - (\nabla \Psi^*) \Psi] \cdot d \mathbf{S} 
$$
You can find this expression in Bransden & Joachain's "Physics of atoms and molecules" (page 69 of 2nd ed). Now, to show that $\Psi$ remain normalized as time goes, we should show that
$$
\oint_S [\Psi^* \nabla \Psi - (\nabla \Psi^*) \Psi] \cdot d \mathbf{S} \to 0
$$
when the surface $S$ embraces, so to speak, all the 3-dimensional space (and if $\int_V |\Psi|^2 dV$ is finite). In chapter one, Griffiths shows this in the one dimensional case (and not in a completely satisfactory way), but what about the more realistic three-dimensional case?
 A: For all practical purpose in physics, the existence of $\int_V |\Psi|^2dV$ comes with $\Psi(x,t)$ going to 0 as $\|x\|$ goes to $+\infty$ whereas the spatial derivatives stay bounded. Thus your last integral goes to 0.
As Valter Moretti pointed out in the comments, it is easy to build a $\Psi$ such that $\int_{\mathbb{R}^3} |\Psi|^2dV$ is finite but $\Psi$ goes to infinity as $\|x\|$ goes to $+\infty$. It is already possible to do that in 1D. See this question on Mathematics Q&A and the handful of answers to get the idea, then blend that with this trick to get smooth functions. Shake and despair at Mathematics… Fortunately, it turns out we physicists have managed to safely ignore such issues without significant backlash.
A: @LucJ.Bourhis's answer supplies the needed information for why $\oint_S \rightarrow 0$ as the size of $S\rightarrow \infty$. Put more concretely, examine the normalization integral in spherical coordinates:
$$N\equiv\int |\Psi|^2 r^2\sin\theta \operatorname{d} r \operatorname{d}\theta \operatorname{d}\phi.$$
In order for $N$ to be finite, $\lim_{r\rightarrow\infty} r^{3+\epsilon} |\Psi|^2 = 0$ for some $\epsilon > 0$ (i.e. the argument of the radial integral has to vanish faster than $1/r$). If you examine the components of the gradient of $|\Psi|^2$ (equal to the argument of $\oint_S$ by the product rule) you'll find that it vanishes like at least the derivative of the upper bound on $|\Psi|^2$. With the upper bound on $|\Psi|^2$ being $\propto r^{-3-\epsilon}$, then the upper bound on $\left|\partial_r |\Psi|^2\right|$ is $\propto r^{-4-\epsilon}$. This question on math.stackexchange is relevant to whether this proof is sound.
Explaining why the angular components are irrelevant is left as an exercise.
A: For a system with Hamiltonian $\hat H$, the time-evolution of an arbitrary state $|\Psi\left(\textbf{r},t\right)\rangle$  is governed by 
$$i\hbar\frac{\partial}{\partial t}|\Psi\left(\textbf{r},t\right)\rangle = \hat H |\Psi\left(\textbf{r},t\right)\rangle \tag{1}\label{SE}.$$
We can take the adjoint of Eq.\eqref{SE} and it reads
$$-i\hbar\frac{\partial}{\partial t}\langle\Psi\left(\textbf{r},t\right)| = \langle\Psi\left(\textbf{r},t\right)| \hat H \tag{2}\label{adjSE}.$$ 
Note that we assume the Hamiltonian to be Hermitian, therefore in the right-hand-side of Eq. \eqref{adjSE}, $\hat H^ \dagger$ was replaced by $\hat H$.
Now we get to the original question namely, the time evolution of the normalization $ \mathcal{N}(t) := \langle\Psi\left(\textbf{r},t\right)|\Psi\left(\textbf{r},t\right)\rangle$ of the wavefunction. Therefore we would like to evaluate the following quantity --
$$\mathcal{\dot N}(t) = \frac{\partial}{\partial t}\langle\Psi\left(\textbf{r},t\right)|\Psi\left(\textbf{r},t\right)\rangle.$$
We proceed as
\begin{align}
\mathcal{\dot N}(t) & = \frac{\partial}{\partial t}\langle\Psi\left(\textbf{r},t\right)|\Psi\left(\textbf{r},t\right)\rangle, \\ \\
 & =  \left(\frac{\partial}{\partial t} \langle\Psi\left(\textbf{r},t\right)|\right)|\Psi\left(\textbf{r},t\right)\rangle + \langle\Psi\left(\textbf{r},t\right)|\left(\frac{\partial}{\partial t} |\Psi\left(\textbf{r},t\right)\rangle\right), \\ \\ 
 & = \frac{1}{-i\hbar}\left(-i\hbar\frac{\partial}{\partial t} \langle\Psi\left(\textbf{r},t\right)|\right)|\Psi\left(\textbf{r},t\right)\rangle + \langle\Psi\left(\textbf{r},t\right)|~\frac{1}{i\hbar}~\left(i\hbar\frac{\partial}{\partial t} |\Psi\left(\textbf{r},t\right)\rangle\right), \\ \\
 & = \frac{1}{-i\hbar}~\left(\langle\Psi\left(\textbf{r},t\right)|\hat H\right)|\Psi\left(\textbf{r},t\right)\rangle + \langle\Psi\left(\textbf{r},t\right)|~\frac{1}{i\hbar}~\left(\hat H|\Psi\left(\textbf{r},t\right)\rangle\right), \mathrm{Using~Eqs. \eqref{SE} \& \eqref{adjSE}} \\ \\ 
 & = \frac{1}{-i\hbar}~\langle\Psi\left(\textbf{r},t\right)|\hat H|\Psi\left(\textbf{r},t\right)\rangle + \frac{1}{i\hbar}~\langle\Psi\left(\textbf{r},t\right)|\hat H|\Psi\left(\textbf{r},t\right)\rangle, \\ \\ 
 & = 0.
\end{align}
Since the time-derivative of the normalization is zero, therefore the normalization remains constant in time. 
Things to note: 


*

*For the normalization to be time-independent, one requires unitary time-evolution of the states $|\Psi\left(\textbf{r},t\right)\rangle$, namely $$|\Psi\left(\textbf{r},t\right)\rangle = \mathcal{U}\left(t,t'\right)|\Psi\left(\textbf{r},t'\right)\rangle.$$ 

*For $\mathcal{U}\left(t,t'\right)$ to be unitary operator of the form $$\mathcal{U}\left(t,t'\right) = \exp\left(-i\hat H (t-t')/\hbar\right),$$ one requires $\hat H$ to be Hermitian. 

*The fundamental reason for the states to stay normalized in a Schrodinger-type time-evolution, is the Hermiticity of the Hamiltonian. Although, there are situations where the states can stay normalized under time-evolution for Hamiltonian which are not Hermitian. However, the time-independence of the normalization is guaranted for a unitary time-evolution, namely for a Hermitian Hamiltonian, under the Schrodinger-type equation. 

