Suppose I'm solving $$\frac{d}{dx}\left(K(x)\frac{du}{dx}\right)=f \text{ in }\Omega,$$ $$u=g \text{ on } \partial\Omega$$where $K(x)$ is smooth  and $$f(x) = \left\{        \begin{array}{ll}            0 & \quad x \neq x_0 \\            1 & \quad x =x_0        \end{array}    \right.$$I tried solving this problem by a finite difference method with f as given, but I noticed some strange behavior in the solution as the number of discrete nodes gets larger.  I'm pretty sure the strange behavior is due to the discontinuity at $x_0$.  So, I propose to replace $f$ by a function $f_\epsilon$ such that $f_\epsilon=0$ everywhere except in the small region $B_\epsilon=\{x\text{ s. t. }|x-x_0|<\epsilon\}$.  In $B_\epsilon$, $f_\epsilon$ is a polynomial such that  $$\begin{array}{ll} f_\epsilon(x_0)=1 & \quad f_\epsilon'(x_0)=0 \\ f_\epsilon(x_0)=0 & \quad f_\epsilon'(x_0-\epsilon)=0 \\f_\epsilon(x_0+\epsilon)=0 & \quad f_\epsilon'(x_0+\epsilon)=0\end{array}$$and $f_\epsilon=0$ throughout the rest of the domain.This gives me a polynomial of degree 5 and ensures $f_\epsilon\in C^1(\Omega)$.  However, I suspect that since finite difference methods seek classical solutions, we probably need more differentiability; that is, $f_\epsilon\in C^2(\Omega)$.  It is not immediately obvious to me what other conditions I can/should impose on $f_\epsilon$ to ensure this extra regularity.  So, I pose the following questions:Does $f_\epsilon$ need to be in $C^2(\Omega)$?  What other conditions can I impose on my polynomial to ensure this extra regularity?Are other (non-polynomial) functions better suited to approximate the discontinuous source term?To observe the true behavior of the original discontinuous source term, should I solve a sequence of regularized problems with $f_\epsilon$ such that $\epsilon\rightarrow 0$?
Regularization of a discontinuous source term in an elliptic pde
pde;finite difference
A more typical regularization would be to use$$f_\epsilon (x) = e^{1-\frac{1}{1-((x-x_0)/\epsilon)^2}},\, {\rm for }\, |x-x_0|<\epsilon$$and zero otherwise.This bump function is $C^\infty$, and you can control its width and height arbitrarily.With regard to your questions, $f$ would typically need to be less regular than your solution $u$, so that's not your problem. And, yes, you would want to solve a sequence of problems, but once you have no grid points inside $|x-x_0| < \epsilon$ besides the one at $x_0$, there's no point in shrinking $\epsilon$ any more. You'd need to refine the mesh at that point. Finally, it's kind of an odd forcing function. Are you sure you didn't mean to use a delta function? 