solving $Lu = f$ numerically where $f$ is discontinuous I am dealing with the second order differential equation
$-y''+p(x)y + q(x)y = r(x)$
on $0 \le x \le 1$ with boundary conditions $y(0)=0$ and $y(1) = 0.$
I am using the finite difference method to obtain a numerical solution for $y$.
However, my $r(x)$ is only continuous almost everywhere since it is defined as:
$r(x) = 0$ when $x \neq 1/2$ and $r(x) = 1$ when $x = 1/2$. Actually, the $1/2$ is arbitrary but it's any point inside the interval $[0,1]$.
Of course, if I run my numerical ODE solver, I will get an answer for $y$ given the discontinuous $r(x)$ that I have. But my question is, am I even allowed to use the finite difference method as usual to solve this problem numerically?
If you're curious as to how the $r(x)$ came about, the ODE above (with $q(x)=0$) is the adjoint problem to my original ode.
can anyone lead me to relevant references?
Thanks!
 A: Yes, you may solve the finite-difference problem, but it is unlikely that you get convergence of any order higher than one. Finite volume methods are advised if coefficients or right-hand side is discontinuous.
Another way may be the following. Let $\mathcal L$ be your differential operator:
$$
\mathcal{L}y \equiv -y''(x) + p(x) y'(x) + q(x) y(x)
$$
and your right hand side has the form
$$
r(x) = a(x) + b(x) \theta(x).
$$
where $\theta(x)$ is the Heaviside theta function. I've moved the discontinuity to $x = 0$, but that should not matter.
Let's search for a solution in form
$$
y(x) = R(x) \theta(x) + S(x).
$$
with $R, S$ being smooth.
Differentiating the solution yields
$$\begin{aligned}
y'(x) &= R'(x) \theta(x) + R(x) \delta(x) + S'(x)\\
y''(x) &= R''(x) \theta(x) + \left[R(x) \delta(x)\right]' + S''(x).\end{aligned}
$$
So
$$
\mathcal{L}y = (\mathcal{L}R)\theta(x) + \mathcal{L}S + p(x)R(x)\delta(x)
-\left[R(x) \delta(x)\right]' = a(x) +b(x) \theta(x).
$$
Since $a(x)\delta(x) = a(0)\delta(x),\; (a(x)\delta(x))' = a(0)\delta'(x)$
$$
(\mathcal{L}R)\theta(x) + \mathcal{L}S + p(0)R(0)\delta(x)
-R(0)\delta'(x) = a(x) +b(x) \theta(x)
$$
the problem reduces to a pair of smooth problems
$$
\mathcal{L}S = a(x)\\
\mathcal{L}R = b(x)\\
R(0) = 0\\
R'(0) = 0.
$$
The condition $R'(0) = 0$ is required for the function $y(x)$ to have a continuous derivative. One can solve the Cauchy problem for $R(x)$ and then solve a boundary problem for $S$ with known boudnary conditions.
Example. Solve
$$
y''(x) + y(x) = \theta(x - 1/2), \quad x \in [0, 1]\\
y(0) = y(1) = 0
$$
Let's make a substitution
$$
y(x) = R(x) \theta\left(x - \frac{1}{2}\right) + S(x).
$$
Solving
$$
R'' + R = 1\\
R(1/2) = R'(1/2) = 0
$$
gives
$$
R = 2 \sin^2 \frac{1-2x}{4}.
$$
Solving 
$$
S'' + S = 0\\
S(0) = y(0) - R(0) \theta(-1/2) = y(0) = 0\\
S(1) = y(1) - R(1) \theta(1/2) = y(1) - R(1) = -2\sin^2 \frac{1}{4}\\
$$
gives
$$
S = -2 \csc (1) \sin^2 \left(\frac{1}{4}\right) \sin x
$$
The final solution is
$$
y = 2 \theta \left(x-\frac{1}{2}\right) \sin ^2\left(\frac{1-2 x}{4}\right)-2 \sin ^2\left(\frac{1}{4}\right) \csc (1) \sin x
$$
