Obtaining both Solutions of Linear, Second-Order ODE using DSolve

Question

Bug introduced in 10.4 or earlier and persisting through 11.3

Reported as CASE:3849226.

While attempting to provide a general answer for question 138321, I obtained a result,

Flatten@DSolve[{D[(1 - R/r) D[Φ[r], r], r] - ((n (n + 1))/r^2) Φ[r] == 0}, Φ[r], r, 
    Assumptions -> n ∈ Integers && n > 0]
{* {Φ[r] -> (r^2 C[1] Hypergeometric2F1[1 - n, 2 + n, 3, r/R])/R^2 + 
             C[2] MeijerG[{{}, {1 - n, 2 + n}}, {{0, 2}, {}}, r/R]} *)

which appears to be incorrect. For instance, plot the two ostensibly independent solutions for n == 2 (and x == r/R).

Plot[Evaluate[{x^2 Hypergeometric2F1[1 - n, 2 + n, 3, x],
    -MeijerG[{{}, {1 - n, 2 + n}}, {{0, 2}, {}}, x]/3} /. n -> 2], {x, 0, 1.2}, 
    PlotRange -> All, Exclusions -> None]

Two problems are evident. First, the two solutions are identical (up to a constant multiplier, here 1/3) for x < 1. Second, neither solution is singular at x == 1, even though it is not difficult to show that one solution should have a logarithmic singularity there. Incidentally,

FunctionExpand[MeijerG[{{}, {1 - n, 2 + n}}, {{0, 2}, {}}, r/R], Abs[r/R] < 1]
(* Hypergeometric2F1[-1 - n, n, 1, 1 - r/R] *)

which I am confident is proportional to r^2 Hypergeometric2F1[1 - n, 2 + n, 3, r/R])/R^2 for all positive integer n. although I have not searched the myriad Hypergeometric Function identities to prove it.

So, my questions are

1. Is this a bug?
2. Is there a DSolve work-around to obtain the second independent solution of the ODE?

Addendum: Correct n == 2 solution

Further insight can be gained by solving the ODE for a specific value of n, which DSolve can handle correctly. For instance, with n == 2,

s2 = With[{n = 2}, FullSimplify@Flatten@DSolveValue[{D[(1 - R/r) D[Φ[r], r], r] -
    ((l (l + 1))/r^2) Φ[r] == 0}, Φ[r], r]]
 (* (3 r^2 (4 r - 3 R) R^5 C[1] - 8 R (-24 r^2 + 6 r R + R^2) C[2] + 
    48 r^2 (4 r - 3 R) C[2] (-Log[r] + Log[r - R]))/(12 R^5) *)

The C are arbitrary constants. For convenience in plotting, set them to C[1] -> - R^-3 and C[2] -> - R^2/(16 Pi).

s2s = FullSimplify[{-Coefficient[s2, C[1]]/R^3, -ReIm@Coefficient[s2, C[2]] R^2/(16 Pi)} 
    /. r -> x R, R > 0 && x > 0]
(* {((3 - 4*x)*x^2)/4, {(1 + 6*(1 - 4*x)*x + 6*(3 - 4*x)*x^2*Re[Log[(-1 + x)/x]])/(24*Pi), 
    Piecewise[{{((3 - 4*x)*x^2)/4, x < 1}}, 0]}} *)

Plot[Evaluate[s2s], {x, 0, 1.2}, PlotRange -> {-.4, .4}, Exclusions -> None]

Comparing this plot with the first plot in the question makes clear that the DSolve solution for arbitrary positive integer n is incomplete. (The first solution is blue, Re of the second solution orange, and Im of the second solution green.) Note that the green curve is not independent of the blue curve and could be eliminated by taking a linear combination of the first and second solutions for 0 < x < 1.

Answer 1

eqn = D[(1 - R/r) D[Φ[r], r], 
     r] - ((n (n + 1))/r^2) Φ[r] == 0;

assume = n ∈ Integers && n > 0;

sol = Assuming[assume, DSolve[eqn, Φ, r][[1]]]

(*  {Φ -> Function[{r}, 
       (r^2*C[1]*Hypergeometric2F1[
                1 - n, 2 + n, 3, r/R])/
           R^2 + C[2]*MeijerG[
             {{}, {1 - n, 2 + n}}, 
             {{0, 2}, {}}, r/R]]}  *)

The solution is valid for Abs[r/R] < 1

(eqn /. sol) // FullSimplify[#, assume && Abs[r/R] < 1] &

(*  True  *)

Or for Abs[r/R] > 1

(eqn /. sol) // FullSimplify[#, assume && Abs[r/R] > 1] &

(*  True  *)

However, the general solution fails for Abs[r/R] == 1

(eqn /. sol) // FullSimplify[#, assume && Abs[r/R] == 1] &

(*  Indeterminate == 0  *)

EDIT: The components are separately solutions

(eqn /. (sol /. C[1] -> 0)) // FullSimplify[#, assume && Abs[r/R] < 1] &

(*  True  *)

(eqn /. (sol /. C[2] -> 0)) // FullSimplify[#, assume && Abs[r/R] < 1] &

(*  True  *)

(eqn /. (sol /. C[1] -> 0)) // FullSimplify[#, assume && Abs[r/R] > 1] &

(*  True  *)

(eqn /. (sol /. C[2] -> 0)) // FullSimplify[#, assume && Abs[r/R] > 1] &

(*  True  *)

The discontinuity comes from the second component when Abs[r/R] == 1

(eqn /. (sol /. C[1] -> 0)) // FullSimplify[#, assume && Abs[r/R] == 1] &

(*  Indeterminate == 0  *)

The first component remains valid for Abs[r/R] == 1

(eqn /. (sol /. C[2] -> 0)) // FullSimplify[#, assume && Abs[r/R] == 1] &

(*  True  *)

Plotting the two independent solutions

Manipulate[
 Plot[Evaluate[
   List @@ (Φ[r] /. sol1) /.
    {C[1] -> c1, C[2] -> c2, 
     r -> x*R, n -> m}],
  {x, -1.1, 1.1},
  Exclusions -> {-1, 1},
  PlotRange -> All,
  AxesLabel -> {Style[ToString[x == r/R, TraditionalForm], 14, Bold], 
    None},
  PlotLegends -> Placed["Expressions", Above]],
 {{c1, 1}, -5, 5, Appearance -> "Labeled"},
 {{c2, 1}, -5, 5, Appearance -> "Labeled"},
 {{m, 2, "n"}, Range[5]}]