(The previous version of this attempted answer was too hasty, resulting in some silly nonsense in part...)
The rigorous interpretations of differential equations or other operator equations with unbounded operators, and/or with rough coefficients, involves issues of self-adjoint extensions of restrictions of symmetric operators, so, in the end, the apparent question of "multiplying" things that aren't meant to be multiplied is illusory.
The case of "singular potential" (an example of "exactly solvable model") such as $H=-\Delta+\delta$ on $\mathbb R$ does appear to be asking to multiply possibly-not-differentiable functions by $\delta$, which is problemmatical. But that "problem" is completely parallel to the "problem" of applying $\Delta$ to functions insufficiently differentiable but that their images under $\Delta$ are outside $L^2(\mathbb R)$, although possibly in some Sobolev space. Thus, the $\delta$ can/should be interpreted as the operator $\delta\otimes\delta$ on test functions given by $(\delta\otimes\delta)u=\delta(u)\cdot \delta$. (So, yes, for continuous $u$ this produces a multiple of $\delta$ by $u(0)$). But since this maps outside $L^2$, it is not a legitimate unbounded operator unless $\delta(u)=0$. So the proper domain of the operator $\delta\otimes\delta$ is test functions vanishing at $0$. This also works for $-\Delta+\delta$ (or, more properly, $-\Delta+\delta\otimes\delta$).
That is, the effect of the "operator" $\delta$ is to impose a boundary condition. At the same time, the Friedrichs extension $T$ of that restriction of $-\Delta+\delta$ does have some "exotic" features: $Tu=f$ if and only if $-\Delta u=f+c\cdot \delta$ for some constant $c$, and the domain of $T$ is $H^2+\mathbb C\cdot u_o$ where $u_o\in H^1$ is a solution of $-\Delta u_o=\delta$$(-\Delta+\lambda_o) u_o=\delta$ for some $\lambda_o$ off the real line. But this formulation skirts the issue of literal pointwise multiplication.
Yes, this seems to define away the issue, but it's not really so. For two point-charges on $\mathbb R$, there are "exotic" eigenfunctions: for example, ($-\Delta+\lambda)u=\delta_a+\delta_b$ has some $L^2$ solutions that meet the implied boundary conditions $u(a)=0=u(b)$, namely, all the fragments of sines and cosines vanishing at the endpoints. So, although $-\Delta$ has no $L^2$ eigenfunctions, the two-point singular potential variant does have, but/and the "multiplication" occurring is not literal pointwise at all.
Yes, solutions are in $H^1$ and the singular potentials are in $H^{-1}$ so application of the latter to the former is legitimate. (My earlier remark about this product as somehow pointwise were both silly and irrelevant, in fact.)
A genuine issue arises in higher dimensions, because Dirac deltas are not in $H^{-1}$ any more. In two dimensions, there are other (less canonical, perhaps) self-adjoint extensions of the restriction of $-\Delta$ to $u$ such that $\delta(u)=0$, and the generic ones of these allow singular potentials in somewhat worse Sobolev spaces. But, still, as the dimension goes up, the problem of maintaining the above interpretation seems to require re-interpreting point-charges as something else, otherwise one does not really quite have a self-adjoint operator on a Hilbert space.