Can a physical wavefunction be non-smooth or even discontinuous? Here's an argument that might support the statement that such a non-smooth wavefunction is not physical:
You cannot add a finite number of smooth functions to get a non-smooth function. By fourier expansion theorem, any function can be expressed as a sum of plane waves (which are smooth with respect to spatial dimensions). Hence you must need an infinite number of smooth functions to get a non-smooth function. Now here's the problem. The different plane waves are momentum eigenfunctions and the non-smooth function is a superposition of these momentum eigenfunctions. If you now try to calculate the expectation value of the momentum, due to the fact that the momentum eigenvalues associated with the momentum eigenfunctions are unbounded from above, the expected value of the momentum could blow up (go to infinity). This is specifically what would make the wavefunction non-physical.
$$\psi(x) = \int \widetilde{\psi}(k) e^{ikx} \mathrm{dk}$$
$$E(k)  = \int |\widetilde{\psi}(k)|^2 k \  \mathrm{dk}$$
But for a general smooth function, how do I know whether the fourier coefficients associated with larger and larger momentum eigenmodes taper off fast enough for the expected value of the momentum to converge?
 A: 
From yuggib's answer: "…we consider physically meaningful only wavefunctions such that there exist a continuous and differentiable representative in its equivalence class. However, also this requirement is not physical."

Not quite. A set of countable pointwise discontinuities may be tolerable, at least at first sight, but there is actually a very good physical reason why wave functions are generally required to be continuous and differentiable, even under typical assumptions on the asymptotic behavior: the average energy, in particular the kinetic energy, may become infinite or undefined. A very simple counter example shows why. 
Consider a (non-relativistic) 1D wavefunction completely confined to the positive semiaxis, with a sharp wavefront at the origin, for instance
$$
\Psi(x) = \theta(x) \psi(x), \;\;\; \psi(0) \neq 0 
$$
where 
$$
\theta(x) = \left\{\begin{array}{c}1, \; \text{for}\;x\ge 0 \\0,\;\text{for}\;x<0\end{array}\right.
$$ 
is the Heaviside step function and $\psi$ is square integrable, with standard  behavior for $x\rightarrow \infty$ (see the related question pointed out by knzhou). The discontinuity poses no problem for normalization, since we can always demand
$$
\int_{-\infty}^\infty{dx\; \Psi^*\Psi} \equiv \int_{-\infty}^\infty{dx\;  \theta(x) \psi^*\psi} = 1,
$$
neither does it interfere much with the average position or momentum, which generally turn out to be finite as long as $\psi$ behaves reasonably at infinity:
$$
\langle \Psi | {\hat x}|\Psi\rangle = \int_{-\infty}^\infty{dx\;\theta(x) x|\psi|^2} < \infty
$$
and
$$
\langle \Psi | {\hat p}|\Psi\rangle = -i\hbar \int_{-\infty}^\infty{dx\; \theta(x) \psi^* \frac{d}{dx}\left(  \theta(x) \psi \right)} =
-i\hbar \int_{-\infty}^\infty{dx\; \left(\theta(x) \psi^* \frac{d \psi}{dx} +\frac{1}{2} \psi^*\psi \frac{d\theta}{dx}\right)}  = \\
= -\frac{i\hbar}{2} \int_{-\infty}^\infty{dx\; \theta(x)\left( \psi^* \frac{d \psi}{dx} - \frac{d \psi^*}{dx}\psi \right)} = \int_0^\infty{dx\; \theta(x){\mathcal J}(x)} < \infty
$$
Note that in the last integral above ${\mathcal J}(x) = -\frac{i\hbar}{2}\left( \psi^* \frac{d \psi}{dx} - \frac{d \psi^*}{dx}\psi \right)$ was identified as the local probability flux, which conveys a nice physical interpretation. 
But regardless of any asymptotic behavior, things are very different for the average kinetic energy, since
$$
\langle \Psi | \frac{{\hat p}^2}{2m}|\Psi\rangle = \int_{-\infty}^\infty{dx\; \Psi^*(x) \left(-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\right) \Psi(x)} = -\frac{\hbar^2}{2m} \int_{-\infty}^\infty{dx\; \theta(x) \psi^*(x) \frac{d^2}{dx^2}\left(  \theta(x) \psi(x) \right)} =
$$
$$
\frac{\hbar^2}{2m} \int_{-\infty}^\infty{dx\; \frac{d}{dx}\left[  \theta(x) \psi^*(x) \right] \frac{d}{dx}\left[  \theta(x) \psi(x) \right]} = \frac{\hbar^2}{2m} \int_{-\infty}^\infty{dx\;  \psi^*\psi \left( \frac{d\theta}{dx}\right)^2} + \text{finite terms} = \frac{\hbar^2}{2m} |\psi(0)|^2\delta(0) + \text{finite terms}
$$
So a sharp wavefront discontinuity, with $\psi(0) \neq 0$, automatically means an uncomfortable infinity in kinetic energy, and in general means bad news for any observable involving 2nd order or higher derivatives.
To answer the question: the average momentum doesn't necessarily blow up in the presence of discontinuities, but the kinetic energy and any averages involving higher derivatives do pose a problem. The asymptotic behavior of the wave function in momentum space must be such that relevant observables involving position derivatives have finite averages. This in turn implies that discontinuities of the type discussed here must be smoothed out.  
A: There is no physical requirement for a wave function to be smooth, or even continuous. At least if we accept the common interpretation that a wavefunction is nothing else than a "probability amplitude". I.e. it represents, when multiplied with its complex conjugate, a probability density. Now a probability density can be discontinuous or even undefined on sets of Lebesgue measure zero. As a matter of fact, the continuous and smooth wavefunction
$$c(x)=e^{-x^2}$$
and the discontinuous wavefunction
\begin{equation}d(x)=\left\{\begin{aligned}&e^{-x^2}&x\neq 0\\&5&x=0 \end{aligned}\right.\end{equation}
give the same probability distribution. In other words, for the purpose of quantum mechanics they are perfectly indistinguishable, and in fact they belong to the same equivalence class of $L^2(\mathbb{R}^d)$ functions. In again different words, the Hilbert space vector $\psi$ on $L^2(\mathbb{R}^d)$ defined by $c(x)$ and $d(x)$ is the same, and therefore it corresponds to the same quantum state.
For the reason above already, continuity is in my opinion meaningless when we deal with wavefunctions. Of course it is always possible (and convenient) to embed, let's say, functions of rapid decrease $\mathscr{S}(\mathbb{R}^d)$ in $L^2(\mathbb{R}^d)$, i.e. considering, by an abuse of notation, $f\in \mathscr{S}(\mathbb{R}^d)\subset L^2(\mathbb{R}^d)$; however we are implicitly considering the equivalence class of almost everywhere equal functions $[f]\in L^2(\mathbb{R}^d)$ to which $f$ belongs, in place of $f$ itself. So let's suppose now that we modifiy OP's assertion to be:

It is unphysical to consider wavefunctions that do not belong to $\mathscr{C}^1_{L^2}(\mathbb{R}^d)=\Bigl\{[f]\in L^2(\mathbb{R}^d), f\in C^1(\mathbb{R}^d)\text{ with }\int_{\mathbb{R}^d}\lvert f(x)\rvert^2 dx<\infty \Bigr\}$.

Now this is a better defined assertion from a mathematical point of view, and it means that we consider physically meaningful only wavefunctions such that there exist a continuous and differentiable representative in its equivalence class. However, also this requirement is not physical. In fact, there are potentials $V(x)$ such that the Hamiltonian $H=-\Delta/2m +V(x)$ is self-adjoint, and $$e^{-iHt}\Bigl[\mathscr{C}^1_{L^2}(\mathbb{R}^d)\Bigr]\nsubseteq \mathscr{C}^1_{L^2}(\mathbb{R}^d)\; .$$
In other words, there are quantum evolutions for which some "smooth" (in the sense just above) wavefunctions become non-smooth as time passes. One relevant concrete example would be the Coulomb potential $V(x)=\pm \frac{1}{\lvert x\rvert}$.
A: Discontinuity in first derivative of wave function implies that the wave function experiences a sudden force that changes its momentum instantly. Thus physically speaking this is not possible as there are no dirac delta potentials. There are potentials very close to dirac delta and thus in the dirac delta approximation the wavefunction will have a first derivative discontinuity.
