Continuity of solutions of ODEs with Dirac delta forcing The problem is : Find $y$, if $y''-4y=\delta(x-a)$, where $\delta$ is the Dirac delta function, and $y$ is bounded as $|x|\rightarrow\infty$.
The solution goes like this:
For $x<a$, $y=Ae^{2x}+Be^{-2x}$, and because of the boundary condition, $B=0$, so $y=Ae^{2x}$. 
For $x>a$, $y=Ce^{2x}+De^{-2x}$, and because of the boundary condition, $C=0$, so $y=De^{-2x}$.
At $x=a$, the solution must be continuous, so $Ae^{2a}=De^{-2a}$. 
Also, the integral from $a^-$ to $a^+$ of the LHS of the equation must be equal to the same integral of the RHS of the equation, which is 1. This gives the other needed equation for $A$ and $D$ from which we can easily find the general solution (I didn't write it as it is irrelevant for the question, but if anyone cares, the general solution is $y=-\dfrac{e^{-2|x-a|}}4$ for $y\neq a$).
My question is: why must $y$ be continuous at $a$? 
An explanation that I have seen is the following:
If $y$ had a simple discontinuity (finite jump) at $a$, then $y'$ would have a $\delta$-function singularity at $a$, and so would $y''$, but there is nothing in the equation to balance such a singularity.
I think that this is an incorrect explanation, as $y''$ would actually have a -2-singularity at $x=a$, while $\delta$-singularity is -1-singularity. This offers an explanation, as the equation doesn't balance out -2-singularities, but is there a more rigorous and simpler proof why $y$ is continuous at $a$?
 A: Indeed, the quoted explanation is incorrect: the singularity of $y''$ would be of higher order (a distribution of order $1$, to be precise), and there is no such thing in the equation. 
The question of "why must the solution be continuous" is always tied to "what do we mean by solution"? For that matter, what does the equation $y''−4y=\delta (x−a)$ mean? What is $y''$ here? The equation cannot be understood in terms of classical derivatives (limit of difference quotient). If $y''$ is the classical 2nd order derivative of $y$, then it is a function, and being added to $-4y$, it is still a function and not a point mass.
One of possible answers is: $y$ is a distribution such that $y''$ is a Radon measure, and the sum of two measures $y''$ and $-4y(t)\,dt$ is the point mass at $a$. But if $y''$ is a Radon measure, then $y'$ is a function of bounded variation on every bounded interval. It follows that $y$ is a Lipschitz function on every bounded interval, and in particular it is continuous. 

Update. One possible answer is: $y$ must be continuous because we say so. It's up to us to decide what we accept as a solution to any problem. For this problem, we decide that $y$ is a solution if it is continuous everywhere, has derivatives $y',y''$ when $x\ne a$, and these derivatives satisfy equation. If such continuous $y$ did not exist, we would then say that the equation has no solution.  
I do not have a satisfactory explanation in terms that you were given. The following is just an attempt to connect the story somehow. Consider the approximate equation where $\delta(x-a)$ is replaced with $D(x-a,h)$. Following the method of doraemonpaul, we find 
$$y=C_1e^{2x}+C_2e^{-2x}+\dfrac{1}{4} \int_a^x \left(e^{2(x-t)}-e^{2(t-x)}\right) D(t-a,h)\,dt$$
As $h\to 0$, $D$ grows large in the sense of pointwise values, but its integral stays equal to $1$. This suggests integration by parts: we throw the integral onto $D(t-a,h)$, turning it into something similar to the Heaviside function, let's say $H(t-a,h)$. 
$$y=C_1e^{2x}+C_2e^{-2x}+\dfrac{1}{2} \int_a^x \left(e^{2(x-t)}+e^{2(t-x)}\right) H(t-a,h)\,dt$$
This is less convenient for practical computation than the original integral, but the advantage is that we have a more reasonable function $H(t-a,h)$ which converges  to a reasonable function $H(t-a)$. In particular, $H$ stays bounded as $h\to 0$, which means that the differences such as $|y(x)-y(x')| $ can be effectively estimated from above:
$$\left|\int_x^{x'} \left(e^{2(x-t)}+e^{2(t-x)}\right) H(t-a,h)\,dt\right|\le M|x-x'|\tag{L}$$ 
where $M$ is the supremum of whatever is under the integral sign. This (L) is the Lipschitz property that I was talking about in the comments, but you can read it simply as continuity. 
A: Just solve it by variation of parameters as usual.
$y_1=e^{2x}$ , $y_2=e^{-2x}$
The Wronskian of these two functions is $\begin{vmatrix}e^{2x}&e^{-2x}\\2e^{-2x}&-2e^{-2x}\\\end{vmatrix}=-4$
$\therefore y=C_1e^{2x}+C_2e^{-2x}-e^{2x}\int_a^x\dfrac{1}{-4}e^{-2x}\delta(x-a)~dx+e^{-2x}\int_a^x\dfrac{1}{-4}e^{2x}\delta(x-a)~dx=C_1e^{2x}+C_2e^{-2x}+\dfrac{e^{2(x-a)}H(x-a)}{4}-\dfrac{e^{2(a-x)}H(x-a)}{4}=C_1e^{2x}+C_2e^{-2x}+\dfrac{(e^{2(x-a)}-e^{2(a-x)})H(x-a)}{4}$
