Numerical solution of ODE with Delta function I want to model a dynamical system of the form
$\frac{\text{d}x}{\text{d}t} = f(x)+nx\delta(\pi(t-0.2)).
$
The problem is that I have a point source which is reoccurring at fixed time steps (say at 0.2,1.2,2.2...). How can I handle this numerically? 
I have tried to figure out solutions and it seems that there are two "easy" approaches:
1) Consider a grid with a fixed step size and define your delta function as 0 on all grid points except the relevant points. However, the influence of the delta function depends strongly on the step size solving this numerically which is why I think this is not correct.
2) Use a Gaussian to represent the Delta-function. This could also work for adaptive step sizes. However, in this case the results depend on the variance of the Gaussian which should be small. Is this a better approach?
If this approach works, I could generate an array consisting of the Gaussians at different time steps which is zero elsewhere and use this with normal ODE-solvers, right?
The third approach is a little bit nasty in the sense that it cannot be used with "normal" ODE solvers. It would be to evaluate f(x) until we approach the point source and take  $y(+\epsilon)=e^ny(-\epsilon)$ (Numerical way to deal with Dirac delta.). This does not depend on the variance or step size. Should I go with this one? I could use normal ODE-solvers until the fixed time step and then just take the exponential. And do you have any references on that?
Thank you very much for your help!
 A: If you have times $a_0<a_1<..$ and coefficients $c_k$ and want to solve
$$
\dot x(t)=f(x(t))+\sum_{k=0}^N c_kδ(t−a_k)
$$
then you can consider the function $u(t)=x(t)-\sum_{k=0}^N c_kH(t−a_k)$ where $H$ is the Heaviside or unit jump/ramp function. Then you get that
$$
\dot u(t)=\dot x(t)-\sum_{k=0}^N c_kδ(t−a_k)=f\left(u(t)+\sum_{k=0}^N c_kH(t−a_k)\right)
$$
which merely has a discontinuous right side. Numerical solvers of the usual sort with adaptive step size control will not be very happy with that situation, but usually can deal with it in a satisfying manner. 

In the modified problem
$$
\dot x(t)=f(x(t))+x(t)\sum_{k=0}^N c_kδ(t−a_k)
$$
you can again try to combine $\dot x(t)-x(t)\sum_{k=0}^N c_kδ(t−a_k)$ into a derivative of some substitution variable. This is achieved in this case with the help of an integrating factor, resulting in setting $x(t)=\exp\left(\sum_{k=0}^N c_kH(t−a_k)\right)u(t)$ where now $u$ satisfies the ODE
$$
\dot u(t)=\exp\left(-\sum_{k=0}^N c_kH(t−a_k)\right)f\left(\exp\left(\sum_{k=0}^N c_kH(t−a_k)\right)u(t)\right).
$$
This has again a piecewise continuous right side, so that the solution $u$ is continuous and piecewise continuously differentiable. Again numerical methods can find this solution, in the post-processing the multiplication with the jumpy exponential introduces the jumps for $x$.
