As proposed by @ZeroTheHero in the long discussion below my question About the traditional explanation of the continuity of the first derivative of a 1D wavefunction, I write down here the rigorous argument popularly attributed to H. Weyl to prove the continuity of the first derivative.
Actually the original argument by Weyl was handled by other people and a relatively more recent source should be Helwig: Differential Operators in Mathematical Physics Adisdon Wesley 1964, Cap 11.
(Unfortunately, I do not have this book, and what I write below is a re-construction of the argument extracted from an Italian theoretical physics textbook.)
Let us consider a ''naive Hamiltonian operator'', where $2m=1$ for shortness,
$$H_0 := -\frac{d^2}{dx^2} + U(x) : C^\infty_0(\mathbb{R}) \to L^2(\mathbb{R}) \tag{0}$$
where $U$ is $C^\infty$ except for a finite number of points $x_k$ where it has just finite jump discontinuities $U(x_k^+)\neq U(x_k^-)$ are finite. (The Schwartz space ${\cal S}(\mathbb{R})$ can be used in place of $C_0^\infty(\mathbb{R})$ achieving the same result.)
According to several results (e.g., by Kato after adding some integrability conditions on $U$ that are irrelevant here), $H_0$ is essentially selfadjoint, i.e., the adjont $H_0^\dagger$ of $H_0$ is selfadjoint
$$H:= (H_0^\dagger)^\dagger = H_0^\dagger\:.$$
Since, in the standard mathematical physics formulation of QM, observables are requested to be properly selfadjoint operators, it is assumed that $H$ is the ''true observable''.
I stress that the domain $D(H)\subset L^2(\mathbb{R})$ is larger than the domain $C_0^\infty(\mathbb{R})$ of $H_0$ and it contains functions which are not smooth. Indeed $H$ is not a differential operator differently than $H_0$.
However $H$ and its domain $D(H)$ are completely determined by $H_0$ and its domain $D(H_0):=C_0^\infty(\mathbb{R})$ through the definition of adjoint operator. In other words, physics is already embodied in $H_0$.
Nevertheless, the existence of a basis of (proper) eigenfunctions (under some standard hypotheses on $U$) is guaranteed for selfadjoint operators as $H$ and not for symmetric operators as the original $H_0$.
So, when dealing with the eigenvector problem we should refer to $H$ and not $H_0$.
From this perspective, the ''correct'' Schroedinger equation is
$$H\psi_E = E\psi_E$$
where $\psi_E \in D(H)$.
It turns out that (there are many theorems leading to this result)
$H$ has the same form as $H_0$ in (0), but the derivatives $\frac{d^2}{dx^2}$ are (second) weak derivatives and they coincide with standard derivatives when $x$ is not a discontinuity point of $U$.
Saying that $g\in L^2(\mathbb{R})$ is the weak derivative (aka distributional derivative) of $\psi \in L^2(\mathbb{R})$ means that
$$\int f(x) g(x) dx = -\int \frac{df }{dx} \psi(x) dx\:, \quad \forall f\in C_0^\infty(\mathbb{R})$$
The Schroedinger equation for $\psi \in L^2(\mathbb{R})$,
$$H\psi = E\psi$$
therefore implies (actually is equivalent to)
$$\int \psi(x)\frac{d^2f}{dx^2} dx = \int (E-U(x)) f(x) \psi(x) dx \:, \quad \forall f\in C_0^\infty(\mathbb{R})\:. \tag{1}$$
Here comes the Weyl result. It states that, if $U$ satisfies the said hypotheses, then
$\quad \quad\quad \quad$ $\psi$ is properly $C^2$ out of the discontinuity points of $U$. $\quad \quad\quad \quad$ [WEYL]
That is a remarkable result as, in principle, $\psi$ is only $L^2$ in this discussion.
The result above has a known pair of fundamental consequences. The latter is the wanted result.
(A) If $\psi$ satisfies (1), then [WEYL] implies that it also satisfies the usual differential Schroedinger equations in the set of points $x\in \mathbb{R}$ where $U$ is continuous.
PROOF.
If $f\in C_0^\infty(\mathbb{R})$ smoothly vanishes on an arbitrary small neighborhood the set of (isolated and finitely many) discontinuity points of $U$, since $\psi$ is $C^2$ where $f$ does not vanish, we can take advantage of the integration by parts obtaioning
$$\int \left(\psi(x)\frac{d^2f}{dx^2} dx - \frac{d^2\psi}{dx^2} f(x) \right) dx = \int \frac{d}{dx}\left(\psi(x)\frac{df}{dx} - \frac{d\psi}{dx} f(x) \right) dx = \psi(b)\frac{df}{dx} - \frac{d\psi}{dx} f(a)=0 $$ where $-a,b>0$ are arbitrary large numbers outside the support of $f$.
Therefore, from (1),
$$\int \psi(x)\frac{d^2f}{dx^2} dx = -\int (U(x)-E) f(x) \psi(x) dx\:,$$
that is, since $\psi$ is $C^2$ where $f$ does not vanish, we can integrate again by parts twice obtaining
$$\int \frac{d^2\psi}{dx^2}f(x) dx = - \int (U(x)-E) f(x) \psi(x) dx\:.$$
Equivalently
$$\int \left(\frac{d^2\psi}{dx^2} - (U(x)-E) \psi(x)\right) f(x) dx= 0\:.$$
The fact that
$-\frac{d^2\psi}{dx^2} + (U(x)-E) \psi(x)$ is continuous outside the discontinuities of $U$ and
arbitriness of $f\in C_0^\infty(\mathbb{R})$ -- exploiting a standard argument of elementary calculus of variations -- implies that, if $x$ is a point where $U$ is continuous:
$$-\frac{d^2\psi}{dx^2} + (U(x)-E) \psi(x)=0$$
(B) Let $\psi$ be as in (A). Then $\psi \in C^1(\mathbb{R})$. In particular it admits continuous first derivative at the discontinuities of $U$.
PROOF. We have to prove that $\psi$ is continuous, differentiable and admits continuous derivative only at the (isolated) points where $U$ is discontinuous, since these facts are already proved by (A) in the remaining points. Let us assume that $x=0$ is a discontinuity of $U$ (where as we know there is a finite jump). From (A) we have that
$$\int_{-\infty}^{0_-} f(x)\frac{d^2\psi}{dx^2} dx + \int_{-\infty}^{0_-} V(x) f(x) \psi(x) dx = E \int_{-\infty}^{0_-}f(x) \psi(x) dx \:, \tag{2}$$
when the support of $f$ is sufficiently narrowed around $0$, so that it does not touch the other discontinuity points of $U$.
Simlarly,
$$\int^{+\infty}_{0_+} f(x)\frac{d^2\psi}{dx^2}+ \int^{+\infty}_{0_+}V(x) f(x) \psi(x) dx = E \int^{+\infty}_{0_+}f(x) \psi(x) dx \:. \tag{3}$$
Summing both sides of the found identities we find
$$\int_{-\infty}^{0_-} f(x)\frac{d^2\psi}{dx^2} dx + \int^{+\infty}_{0_+} f(x)\frac{d^2\psi}{dx^2} dx = \int (E-V(x)) f(x) \psi(x) dx $$
where we used the fact that $0\in \mathbb{R}$ has zero measure and $f(x)\psi(x)(E- U(x))$ is integrable around $x=0$ because $f\psi \in L^1$ and $U$ is bounded on the support of $f$ as it has a just finite jump. Taking (1) into account, the found identity can be arranged to
$$\int_{-\infty}^{0_-} f(x)\frac{d^2\psi}{dx^2} dx + \int^{+\infty}_{0_+} f(x)\frac{d^2\psi}{dx^2} dx = \int \frac{d^2 f}{dx^2} \psi(x) dx $$
Using integration by parts in the left-hand side we find
$$ \left(-\frac{d\psi}{dx}(0_+)+ \frac{d\psi}{dx}(0_-)\right)f(0) + (\psi(0_+)- \psi(0_-))\frac{df}{dx}|_{x=0}+ \int_{-\infty}^{0_-} \psi(x)\frac{d^2f}{dx^2} dx + \int^{+\infty}_{0_+} \psi(x)\frac{d^2f}{dx^2} dx
= \int \frac{d^2 f}{dx^2} \psi(x) dx\:. $$
Namely, since $$ \int_{-\infty}^{0_-} \psi(x)\frac{d^2f}{dx^2} dx + \int^{+\infty}_{0_+} \psi(x)\frac{d^2f}{dx^2} dx
= \int \frac{d^2 f}{dx^2} \psi(x) dx\:,$$
we have
$$ \left(-\frac{d\psi}{dx}(0_+)+ \frac{d\psi}{dx}(0_-)\right)f(0) + (\psi(0_+)- \psi(0_-))\frac{df}{dx}|_{x=0}
= 0 \:.$$
As $f$ is arbitrary (with the said constraints regarding its support), we can choose first $f$ with $f(0)=0$ but $f'(0) \neq 0$ and next $f$ such that $f'(0)=0$ but $f(0) \neq 0$, finding that
$$\psi(0_+)= \psi(0_-)\quad and \quad \frac{d\psi}{dx}(0_+)= \frac{d\psi}{dx}(0_-)$$
as wanted.
