But how can we guarantee that two solutions $\boldsymbol {\psi_1}$ and $\boldsymbol {\psi_2}$ to the time-dependent equation don't have $\boldsymbol {\psi_1(x,0)} = \boldsymbol {\psi_2(x,0)}$. If we can't guarantee this, then how do we know that the solution found by Griffith's method is unique?
I interpret that your question basically asks how do we know that the time-dependent Schrödinger equation with an initial condition has a unique solution. To show this, we use a common method for showing the uniqueness of a solution to a linear differential equation: we show that if two answers are the same at one point, then their difference is constant $0$.
The time-dependent Schrödinger equation is linear; that means that if $\Psi_1$ and $\Psi_2$ are solutions to the time-dependent Schrödinger equation, then so is $\alpha_1 \Psi_1 + \alpha_2 \Psi_2$ for any $\alpha_1, \alpha_2$. In particular, $\Psi=\Psi_1-\Psi_2$ is a solution. So we have
$$
-i\hbar \partial_t \Psi = \hat{H} \Psi,\quad\Psi(0)=0
$$
Now we get
$$
\begin{align*}
\partial_t \langle\Psi, \Psi\rangle
&= \langle\partial_t \Psi, \Psi\rangle + \langle\Psi, \partial_t \Psi\rangle \\
&= \langle \frac{i}{\hbar} \hat{H} \Psi, \Psi \rangle + \langle \Psi, \frac{i}{\hbar} \hat{H} \Psi\rangle \\
= & -\frac{i}{\hbar}\langle\hat{H}\Psi,\Psi\rangle + \frac{i}{\hbar} \langle \Psi, \hat{H} \Psi\rangle \\
&= -\frac{i}{\hbar}\langle\Psi,\hat{H}\Psi\rangle + \frac{i}{\hbar} \langle \Psi, \hat{H} \Psi\rangle \\
&= 0,
\end{align*}
$$
where the next-to-last equation is due to the fact that the Hamiltonian is Hermitian. So now we have a differential equation for the real function $f$ defined by $f(t) \equiv \langle\Psi(t),\Psi(t)\rangle$,
$$
\partial_t f = 0 \, .
$$
Zero derivative means that $f(t)$ is a constant. Therefore,
$$f(t)=f(0)=\left<\Psi(0),\Psi(0)\right>=0 \, .$$
Therefore, $\Psi=0$, so $\Psi_1=\Psi_2$.
You can prove this without looking at any of the specific cases by doing a first-order perturbation of the differential equation that defines the time-evolution operator.
We start with
$$
i \hbar \frac{\partial}{\partial t} U(t, t_0) = H U(t, t_0).
$$
Or, rearranging some terms,
$$
\frac{\partial}{\partial t} U(t, t_0) = - \frac{i}{\hbar} H U(t, t_0).
$$
At a time $t + \delta t$, the above equation tells us that to first order in $\delta t$ we have
$$
U(t + \delta t, t_0) = U(t, t_0) - \frac{i}{\hbar} H U(t, t_0) \delta t.
$$
If $U$ is unitary, then we have
$$
I = U^{\dagger} U (t + \delta t, t_0) = \left[ U^{\dagger}(t, t_0) + \frac{i}{\hbar} U^{\dagger}(t, t_0) H^{\dagger} \delta t\right] \left[ U(t, t_0) - \frac{i}{\hbar} H U(t, t_0) \delta t \right].
$$
Expanding this and keeping only terms to first order in $\delta t$ yields
$$
I = U^{\dagger} U (t, t_0) - \frac{i}{\hbar} U^{\dagger} (t, t_0) H U (t, t_0) \delta t + \frac{i}{\hbar} U^{\dagger}(t, t_0) H^{\dagger} U (t, t_0) \delta t.
$$
Since we demanded that $U$ is unitary for all time, we should also have $U^{\dagger} U (t, t_0) = I.$ Substituting this and canceling $I$ from both sides of the equation yields
$$
0 = \frac{i}{\hbar} U^{\dagger}(t, t_0) \left(H^{\dagger} - H \right) U(t, t_0) \delta t.
$$
We must therefore have $H^{\dagger} = H$. So we have shown that if $U$ is unitary, then $H$ must be Hermitian.
For the other direction, we return to our first order expansion:
$$
U^{\dagger} U (t + \delta t, t_0) = U^{\dagger} U (t, t_0) - \frac{i}{\hbar} U^{\dagger} (t, t_0) H U (t, t_0) \delta t + \frac{i}{\hbar} U^{\dagger}(t, t_0) H^{\dagger} U (t, t_0) \delta t.
$$
If $H$ is Hermitian, then
$$
U^{\dagger} U (t + \delta t, t_0) = U^{\dagger} U(t, t_0).
$$
So $U^{\dagger} U(t, t_0)$ is constant for all $t$. Since $U(t_0, t_0)$ is the identity, we must have $U^{\dagger} U(t, t_0) = I$ for all $t$. We have thus proven that if $H$ is Hermitian, then $U$ must be unitary.
Best Answer
That exponential factor is time evolution factor.
$$|\psi(t)\rangle = e^{-i \hat{H} t/\hbar} |\psi(0)\rangle $$
Such exponential factor comes from unitary property of time evolution operator to maintain its normalization condition. Here, $H$ is Hamiltonian which has eigenvalues as energy.
Energy eigenvalues for given hard wall potential can be obtained from following process:
$$-\frac{\hbar^2}{2m} \frac{d\psi}{dx} = E \psi \; \rightarrow \; \psi \sim \sin kx$$
$$k= \sqrt{\frac{2mE}{\hbar^2}} = \frac{n \pi}{L}$$
$$E_n = \frac{n^2 \pi^2 \hbar}{2mL^2}$$
Above DE has solutions you written, and $k$ is quantized by boundary condition. If you put those energy eigenvalues into time evolution factor, then you can get desired solution.
Linear combination means any wave function in this system can be expanded with energy eigenstates. It's because $\sin nx/L$ has completeness. (Think Fourier series)