plotting - Solving a differential equation involving the square of the derivative

I want to solve the differential equation $$\left( \frac{dr}{d\lambda} \right)^2 = 1-\frac{L^2}{r^2} \left( 1-\frac{1}{r} \right)$$ where $L$ is some parameter. The behavior I am looking for is when $r$ starts "at infinity" at $\lambda=0$, reaches some minimum $r$, and increases out to infinity again. (I'm trying to describe the path that light takes in the presence of a Schwarzschild black hole.)

The issue I'm running into is when I use ParametricNDSolve, I can't tell Mathematica to smoothly transition from negative $dr/d\lambda$ to positive $dr/d\lambda$ after reaching the minimum $r$.

f[r_]:=1-1/r;
soln = ParametricNDSolve[{r'[\[Lambda]]^2 == 1 - L^2/r[\[Lambda]]^2 *f[r[\[Lambda]]], r[0] == 1000}, r, {\[Lambda], 0, 1000}, {L}]

Plot[r[30][\[Lambda]] /. soln, {\[Lambda], 0, 1000}]

In the above code, it plots the graph fine. However if I try to increase the range of $\lambda$ (say to 1200), it breaks down; the situation where I solve for $r'(\lambda)$ first and take the negative square root is identical. I'm not sure how to capture this extra information of transitioning from the negative to positive square root in the differential equation.

Sorry if this is a basic question, I'm rather new to Mathematica (and this site).

Answer

You could differentiate to make a second-order equation without square root problems. Actually, it moves the square trouble to the initial conditions, but that is easier to solve.

foo = D[r'[λ]^2 == 1 - L^2/r[λ]^2*(1 - 1/r[λ]), λ] /. Equal -> Subtract // FactorList
(*  {{-1, 1}, {r[λ], -4}, {r'[λ], 1}, {-3 L^2 + 2 L^2 r[λ] - 2 r[λ]^4 (r^′′)[λ], 1}}  *)

ode = foo[[-1, 1]] == 0     (* pick the right equation by inspection *)

(*  -3 L^2 + 2 L^2 r[λ] - 2 r[λ]^4 (r'')[λ] == 0  *)

icsALL = Solve[{r'[λ]^2 == 1 - L^2/r[λ]^2*(1 - 1/r[λ]), r[0] == 1000} /. λ -> 0,
  {r[0], r'[0]}]
(*
  {{r[0]  -> 1000, 
    r'[0] -> -(Sqrt[1000000000 - 999 L^2]/(10000 Sqrt[10]))}, (* negative radical *)
   {r[0]  -> 1000, 
    r'[0] -> Sqrt[1000000000 - 999 L^2]/(10000 Sqrt[10])}}    (* positive radical *)
*)


ics = {r[0], r'[0]} == ({r[0], r'[0]} /. First@icsALL); (* pick negative solution *)

soln = ParametricNDSolve[
   {ode, ics},
   r, {λ, 0, 2000}, {L}];
Plot[r[30][λ] /. soln, {λ, 0, 2000}]