I am trying to numerically evaluate a Greens function for this equation:$$ \left[\frac{\partial^2}{\partial x^2} + f(x) \right] G(x) = \delta (x-x_0) $$ With Neumann boundary conditions. Here, the function $f(x)$ is known. The difficulty in solving this lies in the numerical approximation of the $\delta$ function. Rather than approximate $\delta \sim \frac{1}{x}$ I attempted to resolve the $\delta$ by integrating both sides in a small interval around $x_0$, but this leaves me with conditions on the derivative of $G$ at $x_0$ that I don't know how to implement in a finite elements method.I would really appreciate a reference for solving such equations (preferably with a finite elements method), or if someone could lend a guiding hand for solving such an equation.
Numerically find Greens function
finite element;finite difference