# DSolve around a jump condition

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.

• Please translate your equations into Mathematica syntax and provide at least your initial attempts at using DSolve etc. Commented Feb 21, 2017 at 22:49
– Karl
Commented Feb 22, 2017 at 2:31
• You need to solve this symbolically? Commented Feb 22, 2017 at 3:55
• What are the other boundary conditions, say, at rmin and rmax? Is l an integer? Commented Feb 22, 2017 at 18:57
• Yes, it needs to be solved symbolically, numerically isn't useful. r_min is zero and r_max is infinity. $l$ is an integer. I've edited the original post for this.
– Karl
Commented Feb 22, 2017 at 20:30

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.

• It's not that I can't find a solution, and the $l$ thing isn't a problem for me, what I want to know is how to code in the boundary conditions that I described.
– Karl
Commented Feb 23, 2017 at 16:16
• @Karl Your ODE falls in the category of "Hybrid", as described in the DSolve documentation. Unfortunately, DSolve returns unevaluated, when I try to solve it as a hybrid equation with WhenEvent. So, I believe that the way I presented is the best available. In general, DSolve should be viewed as "work in progress". Here is what I actually tried. With[{l = 1}, Flatten@DSolve[{D[(1 - R/r) D[Φ[r], r], r] - ((l (l + 1))/r^2) Φ[r] == 0, Φ[0] == 0, WhenEvent[r == 1.1, Φ[r] -> Φ[r] + qq]}, {Φ[r]}, {r, 0, 2}]] Commented Feb 23, 2017 at 16:45