# Need help on working with parametric function

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.

• sol has 2 parameters: u0 and c. Therefore l also should depend on u0 and c. And you try to ContourPlot with two additional variables h and h1, I dont't understand what is your goal. – Alx Oct 31 '19 at 14:09
• @Alx I updated the Q. – all Oct 31 '19 at 17:53
• You write log. Do you really mean Log? – m_goldberg Nov 1 '19 at 3:59
• yes I meant 'Log' – all Nov 1 '19 at 4:14
• @m_goldberg To make the code simple : assume p1 = 10h-l and p2 = 10 h1. – all Nov 1 '19 at 4:38

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}] • Thanks @Alx.I do not see when p1-p2=0 – all Nov 1 '19 at 5:45
• Lines in ContourPlot correspond to p1-p2 == 0 which is the same as 10 h - l1 == 10 h1. In ListPlot3D I plot h vs h1 and l1 as a solution to 10 h - l1 == 10 h1 ->h==h1 - l1/10. – Alx Nov 1 '19 at 7:04
• I was expecting to see intersection of two surfaces !! which indicates where the zero pint of p1- p2. – all Nov 1 '19 at 10:17
• Sorry, don't understand. Let's clarify. p1 and p2 are actually planes in terms of h and h1. You can make this plot: ContourPlot3D[Evaluate[(10 h == 10 h1 - #) & /@ l], {h, 2, 5}, {h1, 2, 5}, {z, -10, 10}] and see that p1 - p2 == 0 generates series of planes, and in cross-section (in h - h1 plane) this gives exactly what I showed in ContourPlot in answer. If it is not what you assumed to be, sorry I have no more ideas. – Alx Nov 1 '19 at 12:49
• Thanks Alx. I got your points. What if z and g function of h and h1 : z=160/h and g=0.0714/h1`. Could you please help me with that? – all Nov 2 '19 at 0:52