ODE's involving distributions or radon measures Well I would like to know an approach (method) to solve "singular" ODE's in a non-formal way.
I seek for a method other (simpler) than   the parameters variation, Laplace transform.
The equation would be 
$$Ly=\delta$$
Suppose that we have solutions of the homogeneous problem.
 A: Consider a linear ODE (the coefficients don't have to be constant, but I assume they are nice enough to not cause trouble): 
$$Ly:=y^{(n)}+a_{n-1}y^{(n-1)}+\dots + a_0 y = \delta \tag1 $$
The delta-function always comes from the highest derivative. For $y^{(n)}$ to contain a delta-function at $0$, the next-to-highest derivative $y^{{n-1}}$ must jump by $1$ at zero. Formally,
$$y^{(n-1)}(0 +)=y^{(n-1)}(0-)+1 \tag2$$
All lower order derivatives must be continuous at $0$, otherwise $y^{(n)}$ will contain a higher order singularity that we don't want. 
$$y^{(k)}(0 +)=y^{(k)}(0-),\quad k=0,\dots,n-2 \tag3$$
The general solution of homogeneous equation $Ly=0$ on the interval $(-\infty,0)$ involves $n$ undetermined coefficients $b_1,\dots,b_n$. The general solution of  $Ly=0$ on the interval $(0,\infty,0)$ also involves $n$ undetermined coefficients, say $c_1,\dots,c_n$. The conditions (2)-(3) are $n$ linear relations between these coefficients. You can use them to eliminate $c_1,\dots,c_n$ from your general solution, leaving it in terms of  $b_1,\dots,b_n$. 

An example:
$$y''+y'-6y=\delta \tag4$$
The general solution of homogeneous equation is $y=b_1e^{2x}+b_2 e^{-3x}$ for $x<0$,   and $y=c_1e^{2x}+c_2 e^{-3x}$ for $x>0$. The relations (2) and (3) take the form 
$$
\begin{split}2c_1-3c_2&=2b_1-3b_2+1 \\
c_1+c_2 &= b_1+b_2
\end{split}
$$
Hence $c_1=b_1+1/5$ and  $c_2=b_2-1/5$. The general solution of (4) is 
$$
y=b_1e^{2x}+b_2 e^{-3x}+\frac15(e^{2x}- e^{-3x})\,H(x)
\tag5$$
where $H$ is the Heaviside function.
