plotting - Help in finding a good guesses for ode

This problem solved in 1 by @bbgodfrey. However, when I changed the parameters, the ode becomes very stiff and I did not know which best guess I need to use for y0. I am looking for :

 ParallelTable[Quiet[p[Z, g, k2, 1.5]], {g, 0.01, 0.03, .02}, {Z, 1000, 4000, 500},{k2, 0.0002, 100, 50}]

Does anyone know a method to estimate a good initial guess? I do not want to just put a random tries every time. I tried -.1 < y0 < -.08,but did not find a solution.

Here is my code:

p[Z0_, g0_, k0_, R0_] := Block[{Z = Z0, 
      g = Rationalize[g0, 0], k2 = Rationalize[k0, 0], ϵ = 10^-4, R = Rationalize[R0, 0]}, 
      ps = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], y[ϵ] == y0, y'[ϵ] == 0, 
        WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, {y, y'[R]}, {r, ϵ], R}, {y0}, 
        Method -> "StiffnessSwitching", WorkingPrecision -> 20];
      sol = FindRoot[Last[ps[y0]], {y0, -1}, Evaluated -> False][[1, 2]];
      tot = 4 π sol NIntegrate[r^2 Exp[-First[ps[sol]][r]], {r, ϵ, R}];
      L = Z/g*tot]
DistributeDefinitions[p];

tab = ParallelTable[Quiet[p[Z, g, 0.0002, 1.5]], {Z, 500, 800, 50}, {g, 0.02, 0.06, .005}];
ListContourPlot[tab, DataRange -> {{.02, .06}, {500, 800}}, PlotLegends -> Automatic, 
    FrameLabel -> {g, Z}, ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]

Answer

Simply stated, a "good" initial guess for y0 in FindRoot is one for which ParametricNDSolve can integrate the ODE in the question from r = ϵ to r = R without errors. For the parameter ranges used in question 210296, -1 or any other negative number roughly in the range {-5, 0} worked well. The present question, however, introduces a much wider set of parameters, some of which accommodate only tiny ranges for the initial guesses. However, if the desired guess can be bracketed in y0, then the bracketed region can systematically be shrunk until a valid guess is obtained.

Specifically, fy in the code block below begins by bifurcating, yc = (yl + yu)/2, the bracketing region, and checks by means of tst = First[ps[yc]]["Domain"][[1, 2]] to determine whether ParametricNDSolve is able to integrate the ODE all the way to R for y0 = yc. If so, fy has found a good y00, which it Throws to p and terminates. Otherwise, it takes advantage of the fact that y'[tst] is negative when yc is more negative than the band of good guesses and is positive when yc is more positive than the band of good guesses in order to reduce the size of the bracketing region by a factor of two, after which the process is iterated until a good guess is found. Typically, less than ten iterations are required to find y00.

p[Z0_, g0_, k0_, R0_] := Block[{Z = Z0, g = Rationalize[g0, 0], k2 = Rationalize[k0, 0], 
        ϵ = 10^-4, R = Rationalize[R0, 0], yl = -8, yu = 0, ps, fy, y00, sol}, 
    ps = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], y[ϵ]] == y0, 
        y'[ϵ] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, {y, y'}, {r, ϵ, R}, 

        {y0}, Method -> "StiffnessSwitching", WorkingPrecision -> 20];
    fy := (Do[yc = (yl + yu)/2; 
        tst = First[ps[yc]]["Domain"][[1, 2]];
        If[tst == R, Throw[yc]]; 
        If[Last[ps[yc]][tst] > 0, yu = yc, yl = yc], {i, 50}]; yc);
    y00 = Which[R == First[ps[yl]]["Domain"][[1, 2]], yl, 
        R == First[ps[yu]]["Domain"][[1, 2]], yu, True, Catch[fy]];
    sol = FindRoot[Last[ps[y0]][R], {y0, y00}, Evaluated -> False][[1, 2]];
    L = Z/g 4 π sol NIntegrate[r^2 Exp[-First[ps[sol]][r]], {r, ϵ, R}]]

Applying p to the parameter set introduced by the present question yields.

DistributeDefinitions[p];
ParallelTable[Quiet[p[Z, g, k2, 1.5]], {Z, 1000, 4000, 500}, {g, 0.01, 0.03, .02}, 
    {k2, 0.0002, 100, 50}]

(* {{{-1.36776*10^12, -14119.8}, {-1.422*10^12, -16956.}}, 
    {{-3.13388*10^12, -33652.5}, {-3.2242*10^12, -40906.5}}, 
    {{-5.62603*10^12, -63057.8}, {-5.79668*10^12, -75910.8}}, 
    {{-8.84774*10^12, -102510.}, {-9.06488*10^12, -122341.}}, 
    {{-1.2798*10^13, -152604.}, {-9.58419*10^10, -179866.}}, 

    {{-1.74782*10^13, -213431.}, {-6.31954*10^10, -248535.}}, 
    {{-2.2886*10^13, -287387.}, {-5.70976*10^10, -327845.}}} *)

The computation requires about ninety seconds on my six-processor computer, with only a small portion of that time used to obtain initial guesses. Note that the iteration over Z, which has most elements, is placed at the leftmost position to maximize the advantage of parallelization. Placing it as in the question increases the runtime to four minutes.

In closing, I would make two points. First, the "guess finder function", fy, fails if {yl, yu} does not bracket at least part of the band of good guesses. If fy does not converge to a good guess, sufficiently decreasing yl should solve the problem. Second, WorkingPrecision may need to be increased for some parameters to improve accuracy.