Short answer: it is technically okay because smooth functions are just subset of the $L^2(\mathbb{R}^3)$.
Long answer: $L^2(\mathbb{R}^3)$ addresses more of the inner product than the differentiability: for normally we want the our eigenfunctions $\{|\psi_j \rangle \}$ are orthogonal to each other, also normalized spacially in $L^2(\mathbb{R}^3)$.
$$
\langle \psi_i |\psi_j\rangle = \int_{\mathbb{R}^3} \overline{\psi_i}(t,x)\psi_j(t,x) \,dx = 0, \quad i\neq j.
$$
Also normalized:
$$
\langle \psi_i |\psi_i\rangle = \int_{\mathbb{R}^3} \overline{\psi_i}(t,x)\psi_i(t,x) \,dx = 1.
$$
Such that
$$
\frac{d}{dt}\langle \psi_i |\psi_i\rangle = 0
$$
If eigenfunctions satisfy above property, then for any normalized wave function $|\psi\rangle$ spanned using inner product has a diagonal form:
$$
|\psi\rangle = \sum^\infty_{i=1} \langle \psi_i |\psi \rangle\, |\psi_i \rangle,
$$
translated to mathematical notation:
$$
\psi(t,x) = \sum^\infty_{i=1} \alpha_i(t) \,\psi_i(t,x),
$$
where
$$
\alpha_i(t) = \int_{\mathbb{R}^3} \overline{\psi_i}(t,x)\psi (t,x)dx,
$$
so that we can have a physical interpretation, namely, the probability amplitude of the system described by wave function $\psi(t,x)$ being at that state $i$ is $\alpha_i$, a.k.a. $|\alpha_i|^2$ is the probability the system at state $i$ at time $t$. All these benefit from the orthonormality of the basis.
The $\psi_i$'s are found by solving the eigenvalue problem:
$$
H\psi_i = E_i \psi_i,
$$
where $E_i$ is the eigenvalue for the operator $H$. Mathematical analogy would be the eigenfunctions of minus Laplacian on a square:
$$
-\Delta u = k^2u,
$$
the eigenfunctions are orthonormal in $L^2$, and orthorgonal in $H^1$.
Differentiability is not the main concern in the quantum world AFAIK.
Firstly, for Bochner space you mentioned, the wave function is normally assumed to be smoothly evolving over time, the obtaining of the wave function for simple models normally follows some Ansatzes, plane wave, decay, etc. The functions we obtain from these Ansatzes are always smooth.
For example, think 1D single particle in a box model:
$$
i\hbar\frac{\partial}{\partial t} \psi(t,x) = H\left|\psi \right> = -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} \psi(t,x) \tag{$\star$}
$$
there is a particle at some state $E$ in $x\in [-1,0]$, then at some time $t_0$ the box length is doubled, $x\in [-1,1]$. What we do is to expand the wave function $\psi(t_0,x)$ at $t_0$ in the old box using the new eigenfunctions $\phi_i$ in the new box $[-1,1]$:
$$\psi(t_0,x) = \sum_{i=1}^{\infty} \langle\phi_i(x)|\psi(t_0,x)\rangle \phi_i(x),\tag{1}$$
where the new eigenfunctions $\phi_i$ solves the eigenvalue problem for $(\star)$ in the new box
$$
-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\phi_i(t,x) = E_i \phi_i(t,x) , \quad \text{ for }x \in[-1,1].
$$
Then
$$
\psi(t_0+\Delta t,x) = \sum_{i=1}^{\infty} \langle\phi_i(x)|\psi(t_0,x)\rangle \phi_i(x)e^{ -iE_n \Delta t/\hbar } \tag{2}
$$
The question is: Does the wave function evolve smoothly? The answer is yes, well, methodologically, there is no non-smoothness existing. For (1) can be expanded in old eigenfunctions or in new eigenfunctions. An exercise for you is to check what happened when $\Delta t\to 0$ in (2).
Second, the differentiability in space is normally infinite within the domain of interest, and continuous up to the boundary of the domain of interest, zero outside. Just think of that particle in the box model (infinite potential wall): the eigenfunction is set to be zero on boundaries $x=1,-1$.