DSolve around a jump condition

Question

I am trying to include a jump condition inside my differential equation, but I'm not exactly sure how to input the boundary conditions. My DE is the following

$\partial_r(f \partial_r \Phi(r))-\frac{l(l+1)}{r^2}\Phi(r)=0$

Which I want to satisfy the condition:

$\Phi(r)=\Phi_>(r)\Theta(r-r_0)+\Phi_<(r)\Theta(r_0-r)$

($\Theta(r)$ being the heaviside or step function)

and jump conditions:

$\Phi_>(r_0)-\Phi_<(r_0)=4\pi q\sqrt{\frac{2l+1}{4\pi}}$

$\Phi'_>(r_0)-\Phi'_<(r_0)=0$

where $r_0$ is an arbitrary point, and $f=1-\frac{R}{r}$ with R being a separate arbitrary point less than $r_0$. $\Phi_<$ refers to the value of the function when $r$ is less than $r_0$. $\Phi_>$ refers to the value of the function when r is greater than $r_0$. All other values should be included. q in this case is just the charge of a particle, l is an integer. The other boundary conditions are such that I can't have the solution blowing up at $r_0$, which I can deal with after I have a solution. I also don't want the solution to blow up at infinity, but that is also something that can be dealt with after the solution.

To be clear, I am working to solve for $\Phi(r)$ in terms of $r$, and want to know how to implement the boundary conditions described by the jump conditions.

f[r_] = 1 - R/r;
DSolve[D[f[r] D[Φ[r], r], r] - ((l (l + 1))/r^2) Φ[r] == 0, Φ[r], r]

The code works fine as is, and gives a hypergeometric function which is expected, but the function isn't of the right form since I haven't included any boundary conditions, which I have no idea how to implement.

Answer 1

I am unable to solve this system for arbitrary positive integer l, perhaps because there is a bug in DSolve. (See question 138440.) However, it can be solved for any particular l. For instance, with l == 1,

s = With[{l = 1}, Flatten@DSolveValue[{D[(1 - R/r) D[Φ[r], r], r] - 
    ((l (l + 1))/r^2) Φ[r] == 0}, {Φ[r], Φ'[r]}, r]]
(* {r^2 C[1] + (C[2] (2 r R + R^2 - 2 r^2 Log[r] + 2 r^2 Log[r - R]))/(2 R^3), 
 2 r C[1] + (C[2] (-2 r + (2 r^2)/(r - R) + 2 R - 4 r Log[r] + 4 r Log[r - R]))/(2 R^3)} *)

This solution is valid for r < r0 and for r > r0 but with different values of C. So, perform the replacement {C[1] -> c1m, C[2] -> c2m} for r < r0, and {C[1] -> c1p, C[2] -> c2p} for r > r0. Also, note that, since r0 > R, c2m must be zero to eliminate the term singular at R. Next, solve for two of the three remaining constants, using the jump conditions in the question.

Flatten@Solve[Thread[(s /. {C[1] -> c1m, C[2] -> 0}) + {qq, 0} == 
    (s /. {C[1] -> c1p, C[2] -> c2p})] /. r -> r0, {c1p, c2p}] // Simplify
(* {c1p -> (R (-qq R + 2 qq r0 + c1m R^2 r0) + 2 qq (R - r0) r0 Log[r0] + 
    2 qq r0 (-R + r0) Log[-R + r0])/(R^3 r0), c2p -> 2 qq (R - r0)} *)

where qq == 4 Pi q Sqrt[(2l+1)/(4Pi)] has been introduced for compactness.

Using other values of l in the computation of s yields analogous but lengthier expressions.