Need help on working with parametric function

Question

Here is an ODE solved by ParametricNDSolveValue solver.I got the solution. Now, I need to solve these two equations:

p1 = 10 h - l;
p2 = 10 h1; 

where h and h1 go from 2 to 5. Left to me. l is obtained from solving s1, where

l = 4 π j NIntegrate[r^2 Exp[-First[s1[sol]][r]], {r, 0, R}];

So, I need to get p1 for different values of l. To do that I need to vary z an g; let's say {z, 100, 400, 50} and {g, 0.023, 0.134}. Then for these values, there is a corresponding l to pass to p1. After getting p1 and p2, I want to know where p1 - p2 == 0 by means of inspecting a 3D contour plot. See the following link 1.

Z = 800; 
g = Rationalize[0.0238, 0];
k2 = Rationalize[0.000194519, 0];
ϵ = 10^-4; 
R = Rationalize[1.5472, 0];
j = .00001271;

s1 = 
  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 -> 30];

sol = FindRoot[Last[s1[y0]][R], {y0, -11.25, -11.}, 
   Evaluated -> False][[1, 2]]

Plot[First[s1[sol]][r], {r, ϵ, R}, 
  AxesLabel -> {r, y}, 
  ImageSize -> Large, 
  LabelStyle -> {Black, Bold, Medium}]

Plot[
  Evaluate @ Table[TooltipFirst[s1[sol]][r], y0] 
  {z, 100, 400, 50}, {g, 0.023, 0.134}, {r, 0, R}] 

l = 4 π j NIntegrate[r^2 Exp[-First[s1[sol]][r]], {r, 0, R}];
p1 = 10 h - l1 ;
p2 = 10 h1 ; 

ContourPlot3D[
  {p1 == 0, p2 == 0}, {h, -2, 2}, {h1, -2, 2}, {l, -2, 2},  
  MeshFunctions -> {Function[{h, h1, l, f}, p1 - p2]},  
  MeshStyle -> {{Thick, Blue}}, 
  Mesh -> {{0}},  
  ContourStyle ->  Directive[Orange, Opacity[0.5], Specularity[White, 30]]]

My issues that l depends on the s1 and therefore on z, g, etc. I believe I need to use Do or Table, but I do not know how? Please help.

Answer 1

I'm not sure I understood your task, but my thought on this.

z = 800; g = Rationalize[0.0238, 0];
k2 = Rationalize[0.000194519, 0];
ϵ = 10^-4; R = Rationalize[1.5472, 0];
j = .00001271;

s1 = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], 
    y[ϵ] == y0, y'[ϵ] == 0, 
    WhenEvent[r == 1, y'[r] -> y'[r] + z1 g1]}, {y, 
    y'}, {r, ϵ, R}, {y0, z1, g1}, 
   Method -> "StiffnessSwitching", WorkingPrecision -> 30];

So, I added two more parameters z1 and g1 (renamed to not to confuse with some initial constants z and g).

Then I define helper function for use inside FindRoot, it accepts numerical arguments only:

s2[y0_?NumericQ, z1_?NumericQ, g1_?NumericQ] := Last[s1[y0, z1, g1]]

Now we can compute y0 values for different z1 and g1:

tab = Flatten[
  Quiet@Table[{y0 /. FindRoot[s2[y0, z, g][R], {y0, -10}], z, g},
  {z, 100, 400, 50}, {g, {0.023, 0.134}}], 1]

and corresponding l values:

l = 4 π j NIntegrate[r^2 Exp[-First[s1[##]][r]], {r, 0, R}] & @@@ tab

What to do next, how to visualize? One can make ContourPlot:

ContourPlot[
 Evaluate@Table[10 h - l1 == 10 h1, {l1, l}], {h, 2, 5}, {h1, 2, 5},
 FrameLabel -> {h, h1}]

Or ListPlot3D:

ListPlot3D[
 Flatten[Table[{h1, l1, h1 + l1/10}, {h1, 2, 5, 0.1}, {l1, l}], 1], 
 AxesLabel -> {h1, l1, h}]

EDIT

To answer additional request: z = 160/h, g = 0.0714/h1. Make corresponding changes:

s1 = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], 
        y[ϵ] == y0, y'[ϵ] == 0, 
        WhenEvent[r == 1, y'[r] -> y'[r] + 160/h 0.0714/h1]}, {y, 
        y'}, {r, ϵ, R}, {y0, h, h1}, 
       Method -> "StiffnessSwitching", WorkingPrecision -> 30];

s2[y0_?NumericQ, h_?NumericQ, h1_?NumericQ] := Last[s1[y0, h, h1]]
ff[h_?NumericQ, h1_?NumericQ] := y0 /. FindRoot[s2[y0, h, h1][R], {y0, -10}]
ll[h_?NumericQ, h1_?NumericQ] := 4 π j NIntegrate[r^2 Exp[-First[s1[ff[h, h1], h, h1]][r]],{r, 0, R}]

Now we are ready to plot contour of 10 h - ll[h, h1] - 10 h1 == 0:

ListContourPlot[
 Flatten[Table[{h, h1, 10 h - ll[h, h1] - 10 h1}, {h, 2, 5, 0.1}, {h1,
     2, 5, 0.1}], 1], Contours -> {0}, ContourShading -> None, 
 ContourStyle -> Thick, FrameLabel -> {h, h1}, BaseStyle -> 12]

FINAL EDIT

k2 = Rationalize[0.000194519, 0];
ϵ = 10^-4; R = Rationalize[1.5472, 0];
j = .00001271;

s1 = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], 
    y[ϵ] == y0, y'[ϵ] == 0, 
    WhenEvent[r == 1, y'[r] -> y'[r] + 160/h 0.0714/h1]}, {y, 
    y'}, {r, ϵ, R}, {y0, h, h1}, 
   Method -> "StiffnessSwitching", WorkingPrecision -> 30];

s2[y0_?NumericQ, h_?NumericQ, h1_?NumericQ] := Last[s1[y0, h, h1]]
ff[h_?NumericQ, h1_?NumericQ] := 
 y0 /. FindRoot[s2[y0, h, h1][R], {y0, -10}]
ll[h_?NumericQ, h1_?NumericQ] := 
 4 π j NIntegrate[
   r^2 Exp[-First[s1[ff[h, h1], h, h1]][r]], {r, 0, R}]

To speed-up computation (FindRoot and NIntegrate are very expensive for continuous plotting) I suggest using Interpolation of ll:

llinterp = 
 Interpolation[
  Flatten[Table[{{h, h1}, ll[h, h1]}, {h, 2, 5, 0.1}, {h1, 2, 5, 
     0.1}], 1]]

Now we can define p1 and p2 via llinterp (as of last comment of OP):

p1 = 10 h - llinterp[h, h1] + u
p2 = -u + 10 h1 llinterp[h, h1]

and see what are ContourPlot3Ds:

ContourPlot3D[p1, {h, 2, 5}, {h1, 2, 5}, {u, 2, 10}, Contours -> 5, 
PlotLegends -> Automatic, AxesLabel -> {h, h1, u}]

ContourPlot3D[p2, {h, 2, 5}, {h1, 2, 5}, {u, 2, 10}, Contours -> 5, 
PlotLegends -> Automatic, AxesLabel -> {h, h1, u}]

Colors in lenends correspond to values of p1, p2.

Apparently, with given ranges of variables (h, h1, u) we can not fulfil conditions p1 = 0 and p2 = 0. Moreover, these conditions are conflicting: to decreese p1 one has to decrease h, but for p2 situation is opposite. So, OP should think about conditions and/or variables' ranges.

We can also plot intersection (p1 == p2):

ContourPlot3D[p1 == p2, {h, 2, 5}, {h1, 2, 5}, {u, 2, 10}, 
AxesLabel -> {h, h1, u}]