So I see your whole question as this:
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?
This is a really basic property of vectors and matrices; if you have two $n\times 1$ "column vectors" $u$, $v$ in $\mathbb C^n$ which are eigenvectors of an $n\times n$ Hermitian matrix $H = H^\dagger$ then either they are orthogonal in the sense that $u^\dagger v = v^\dagger u = 0$ or else they have the same eigenvalue.
Why? Well it's simple interplay of the pieces involved. Hermitian matrices have to have real eigenvalues because $u^\dagger ~ H ~ u = \lambda_u ~ u^\dagger ~ u$ needs to be its own complex conjugate -- because $(A B)^\dagger = B^\dagger A^\dagger$ and the fact that the $1\times 1$ adjoint operation is a complex conjugation you find $(u^\dagger ~ H ~ u)^* = \lambda_u^* ~ u^\dagger ~ u$ is, by the very definition of Hermitian ($H = H^\dagger$), also $ u^\dagger ~H^\dagger~u = u^\dagger ~ H ~ u = \lambda_u u$, hence we either trivially have $u^\dagger ~ u = 0$ and $u = 0$, or else $\lambda_u = \lambda_u^*$ and hence $\lambda_u$ is a real number.
Now look at $u^\dagger \cdot H \cdot v = \lambda_v u^\dagger v$ but its complex conjugate must be $v^\dagger \cdot H \cdot u = \lambda_u v^\dagger u$; taking a second complex conjugate we find that $\lambda_v u^\dagger v = \lambda_u u^\dagger v$. We only have two options here: either $\lambda_u = \lambda v$ or else $u^\dagger v = 0$.
Essentially the problem is: matrices always have the same "left" and "right" eigenvalues, since the expression $\det (A - \lambda I) = 0$ that you use to find them uses a determinant, and $\det A = \det A^T$ because determinants don't care whether you iterate on minors by-row or by-column. When you say that a matrix is its own transpose, as Hermitian ones basically are, then you're saying that also the left- and right-eigenvectors are, respectively, transposes of each other. And this fact forces those eigenspaces to become orthogonal to each other.
Now I claim that the above points do not depend on the discrete indices of the row vector or column vector, but only on the relationship of the objects, and hence when we go from a discrete-indexed vector to a continuous-indexed wavefunction, we get the same conclusion. Let's follow it through. Let $\mathcal A$ map wavefunctions to wavefunctions, so that we have$$\langle \phi | \hat A | \psi \rangle = \int_{\mathbb R^n} d^n r ~ \phi^*(\vec r) ~ \mathcal A[\psi](\vec r).$$ Furthermore let's define the shorthand notation that $$ \langle \phi | \psi \rangle = \langle \phi | 1 | \psi \rangle = \int_{\mathbb R^n} d^n r ~ \phi^*(\vec r) ~ \psi(\vec r).$$We assert that $\mathcal A$ is Hermitian, which means that: $$ \langle \phi | \hat A | \psi \rangle = (\langle \psi | \hat A | \phi \rangle)^* = \left[ \int d^n r ~ \psi^*(\vec r) ~ \mathcal A[\phi](\vec r) \right]^*.$$ It follows that any eigenvalue $\hat A | a \rangle = a | a \rangle$ is real, because we find similarly to the above that $$(\langle a | \hat A | a \rangle)^* = a^* \langle a | a \rangle = \langle a | \hat A | a \rangle = a \langle a | a \rangle, ~~\text{ so }a = a^*.$$ But similarly we have that if $\hat A | a \rangle = a | a \rangle$ and $\hat A | b \rangle = b | b \rangle$ then $$\langle a | \hat A | b \rangle = b \langle a | b \rangle = (\langle a | \hat A | b \rangle)^{**} = (\langle b | \hat A | a \rangle)^* = (a \langle b | a \rangle)^* = a^* \langle a | b \rangle.$$This again forces us to conclude that either $b = a$ or else $\langle a | b \rangle = 0$.
Now consider the very special 1D case where $$\mathcal A[\psi](x) = i~\psi'(x)$$. This imaginary-single-derivative operator is Hermitian, I claim. To see why, let's look at the definition; it's Hermitian if and only if $$\int_{-\infty}^\infty dx ~ \phi^*(x) ~ i~\psi'(x) = \left[ \int_{-\infty}^\infty dx ~ \psi^*(x) ~ i ~ \phi'(x) \right]^*.$$ Is this true? Yes, we can integrate the left hand side by parts, raising $\psi'$ while we lower $\phi^*$, to find: $$\int_{-\infty}^\infty dx ~ \phi^*(x) ~ i~\psi'(x) = \left[i ~ \phi^*(x) ~\psi(x)\right]_{-\infty}^\infty - \int_{-\infty}^\infty dx ~ [\phi^*]'(x) ~ i~\psi(x), $$ the first term of which can be discarded as 0 since both wavefunctions tend to 0 at infinity. This leaves only $$ \int_{-\infty}^\infty dx ~ [\phi'(x)]^* ~ (-i) ~\psi(x) = \int_{-\infty}^\infty dx ~ [\phi'(x) ~ i ~\psi^*(x)]^* = \left[\int_{-\infty}^\infty dx ~ \psi^*(x)~ i ~ \phi'(x)\right],$$which is what we set out to prove.
It follows that $- \frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2}$, the composition of two such operators, which has been multiplied by a constant, is also Hermitian, as is any constant real function -- thus so is $- \frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x)$.
Since it's Hermitian, when we solve $E \psi(x) = - \frac{\hbar^2}{2m} \frac{\partial^2\psi}{\partial x^2} + V(x)~\psi(x)$ we find that any two distinct solutions $\phi, \psi$ either have the same energy or else $\int dx ~ \phi^*(x) ~ \psi(x) = 0$, precluding them from ever being the same function.
Furthermore, if we ever find the complete spectrum $\{\phi_i\}$, this gives us an easy way to find the components:$$\psi(x) = \sum_i c_i \phi_i(x) \text{ where } c_i = \int dx~\phi_i^*(x) \psi(x).$$
