Your solution in the second case returned complex numbers when evaluated. When I switched the method to StiffnessSwitching, the complex number went away. Note that I changed to ParametricNDSolveValue, because that is what I usually work with.
pr2 = ParametricNDSolveValue[{(1 + x)^5 D[(r[x])/(1 + x)^4, x] ==
l024 (r[x])^(1/2), r[zp] == fot}, r, {x, 0, 10^8}, {l024},
WorkingPrecision -> 75, Method -> "StiffnessSwitching"];
ab11 = ContourPlot[((pr2[l024][x])/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2}, PlotLegends -> Automatic, PlotRange -> All]
ab22 = ContourPlot[((pr2[l024][x])/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2},
PlotLegends ->
BarLegend[Automatic, LegendMarkerSize -> 180,
LegendFunction -> "Frame", LegendMargins -> 5,
LegendLabel -> "\!\(\*SubscriptBox[\(z\), \(Lss\)]\)"],
Frame -> True,
FrameLabel -> {{"\!\(\*SubscriptBox[\(z\), \(Lss\)]\)",
""}, {"\!\(\*SubscriptBox[\(\[Lambda]\), \(0\)]\)", ""}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 14},
Contours -> {5}, ContourStyle -> Directive[Thick, Black],
ContourShading -> None, PlotRange -> Full];
Show[ab11, ab22]
ab22
Additional Analysis and Scaling
In the following, I will do some basic analysis and scaling of the differential equation. I will the subscript $d$ to denote the variable/parameter has dimensions. Here is OP's initial equation:
$${\left( {{x_d} + 1} \right)^5}\frac{\partial }{{\partial {x_d}}}\frac{{{r_d}\left( {{x_d}} \right)}}{{{{\left( {{x_d} + 1} \right)}^4}}} = {\lambda _d}\sqrt {{r_d}\left( {{x_d}} \right)} ;{x_d} \geq 0$$
We can use Mathematic to evaluate and simplify the equation to obtain:
$$\frac{{\partial {r_d}\left( {{x_d}} \right)}}{{\partial {x_d}}} = \frac{{{\lambda _d}\sqrt {{r_d}\left( {{x_d}} \right)} + 4{r_d}\left( {{x_d}} \right)}}{{\left( {{x_d} + 1} \right)}}$$
We can define dimensionless variables and parameters like so:
$$x = \frac{{{x_d}}}{{{z_p}}};r = \frac{{{r_d}}}{{{f_{ot}}}};\lambda = \frac{{{\lambda _d}}}{{\sqrt {{f_{ot}}} }}$$
Now, we can create a non-dimensionalized equation like so:
$$\frac{{dr}}{{dx}} = \frac{{4r + \lambda \sqrt r }}{{\frac{1}{{{z_p}}} + x}}$$
We know that at $r(x=1)=1$, which implies that the right hand side of the equation is real and positive. Beyond $x=1$, $r$ is a monotonically increasing function. If we look backwards from $x=1$, then $r$ should be monotonically decreasing. A singularity occurs at $x=-\frac{-1}{z_d}$, but we are always above that point since $x \geq 0$. Examining the equation in simplified non-dimensional form, it is difficult to see how $r$ could turn complex since the right hand side should be positive.
Here is an example workflow using the non-dimensionalized form. I increased the MaxRecursions in the plot to eliminate the small spikes. Also, I imported the NDSolveUtilities package to look at the timesteps taken by the solver.
Needs["DifferentialEquations`NDSolveUtilities`"];
eq = r'[x] == (4 r[x] + λ Sqrt[r[x]])/(1/zp + x);
pr3 = ParametricNDSolveValue[{eq, r[1] == 1},
r, {x, 0, 2}, {λ}, WorkingPrecision -> 75,
Method -> "StiffnessSwitching"];
ab111 = ContourPlot[((pr3[λd/Sqrt[fot]][xd/zp])/(lu/fot))^(1/
4), {λd, 0, 1.2*10^-22}, {xd, 0, 2},
PlotLegends -> Automatic, PlotRange -> All]
ab222 = ContourPlot[((pr3[λd/Sqrt[fot]][xd/zp])/(lu/fot))^(1/
4), {λd, 0, 1.2*10^-22}, {xd, 0, 2}, MaxRecursion -> 4,
PlotLegends ->
BarLegend[Automatic, LegendMarkerSize -> 180,
LegendFunction -> "Frame", LegendMargins -> 5,
LegendLabel -> "\!\(\*SubscriptBox[\(z\), \(Lss\)]\)"],
Frame -> True,
FrameLabel -> {{"\!\(\*SubscriptBox[\(z\), \(Lss\)]\)",
""}, {"\!\(\*SubscriptBox[\(λ\), \(0\)]\)", ""}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 14},
Contours -> {5}, ContourStyle -> Directive[Thick, Green],
ContourShading -> None, PlotRange -> All];
Show[ab111, ab222]
ab222
StepDataPlot[pr3[(1.2*10^-22)/(2 Sqrt[fot])]]
With the StiffnessSwitching method activated, we see a nice smooth transition to the timestep. The following plots show the timestep control for 4 cases that I ran.
Setting the AccuracyGoal only looks like a coarse description of when the StiffnessSwitching is turned on. The WorkingPrecision only setting appear to give up on adjusting the timestep when the solution moves away from the initial boundary condition.
Let's check the assumptions of the previous analysis that said r was monotonically increasing and positive by plotting r vs x and r(0) vs $lambda_d$ with the following code:
Plot[((pr3[0.6*10^-22/Sqrt[fot]][xd/zp])/(lu/fot))^(1/4), {xd, 0,
2 zp}]
Plot[((pr3[λd/Sqrt[fot]][0/zp])/(lu/fot))^(1/4), {λd,
0, 1.2*10^-22}]
The results seem to be consistent with our previous statements.
Finally, let's compare the "ab2" plots of AccuracyGoal Only (red), WorkingPrecision++StiffnessSwitching (green), and WorkingPrecision+AccuracyGoal+StiffnessSwitching (blue).
Show[ab2, ab222, ab2222]
The blue curve took the longest, but had the most control and probably the most accurate. One needs to determine if the extra cost it worth it.
ab22refer toprversuspr2? Is it a potential typo? $\endgroup$StiffnessSwitchingand it seemed to help. $\endgroup$