How to effectively guess the initial values for a Poincaré section?

I have the following code, which plots the Poincaré section for the trajectory of particles around a symmetric black hole:

kr = 100; kθ = 25; rh = 2; rc = 3.2; θc = 0; κ = \
0.25;

f[r_] := 2 κ (r - rh)

H[r_, pr_, θ_, pθ_] := -Sqrt[1 - f[r]] pr + Sqrt[
pr^2 + pθ^2/r^2] + 1/2 kr (r - rc)^2 +
1/2 kθ (rh θ - rc θc)^2
H0[{r_, pr_, θ_, pθ_}] := -Sqrt[1 - f[r]] pr + Sqrt[
pr^2 + pθ^2/r^2] + 1/2 kr (r - rc)^2 +
1/2 kθ (rh θ - rc θc)^2

{a1, b1, c1,
d1} = {D[H[r[t], pr[t], θ[t], pθ[t]], pr[t]],
D[H[r[t], pr[t], θ[t], pθ[t]],
pθ[t]], -D[H[r[t], pr[t], θ[t], pθ[t]],
r[t]], -D[
H[r[t], pr[t], θ[t], pθ[t]], θ[
t]]}; abcd = {r'[t] == a1, pr'[t] == c1, θ'[t] == b1,
pθ'[t] == d1};

ClearAll[solution];
solution[R_, PR_] :=
Reap@NDSolve[
Flatten[{abcd, {r[0] == R, pr[0] == PR, θ[0] == 0,
pθ[0] == initial[R, PR]},
WhenEvent[θ[t] == 0 && θ'[t] > 0,
Sow[{r[t], pr[t]}]]}], {}, {t, 0.0, T}, MaxSteps -> Infinity]

T = 4000.0;

energy = 75.0;

constraint = (pθ /.
Solve[energy == H[r, pr, 0, pθ], pθ]) // Last;

ClearAll[initial];
initial[r_?NumericQ, pr_?NumericQ] := Evaluate[constraint]

tab = Table[solution[R, 0.0] // Last // First, {R, 2.72, 3.75, 0.01}];

ListPlot[tab, PlotRange -> All, PlotTheme -> "Detailed",
PlotStyle -> PointSize[.003], AspectRatio -> 1/GoldenRatio,
ImageSize -> Large, FrameLabel -> {r, Subscript[p, r]},
LabelStyle -> Directive[Black, 10], GridLines -> None]

which results in the desired plot:

In this case, I was reproducing results from a published paper, so I was easily able to input the desired range of R as

{R, 2.72, 3.75, 0.01}, since it was in the paper. But suppose I was making the Poincaré section for the first time, and let's say I guess the range to be {R, 2.01, 4, 0.01}:

tab = Table[solution[R, 0.0] // Last // First, {R, 2.01, 4, 0.01}]

which gives me warning messages and some missing values.

also, it will not give me any plot whatsoever.

So, is there any way using Mathematica by which I can get the desired range?

Edit: Applying the tricks provided by @E. Chan-López and @Daniel Huber, also taking into consideration the useful comment by @Alex Trounev, I tried to scan for desired range using:

tab = Quiet@
Table[If[Flatten[solution[R, 0.0]] === {}, Nothing,
solution[R, 0.0] // Last // First], {R, 2.01, 4, 0.1}]

tmp = Flatten[tab, 1]

MinMax /@ {tmp[[All, 1]], tmp[[All, 2]]}

which gives me the range as:

{{0.594091, 3.90931}, {-26.3422, 868.724}}

It's still not quite correct. I even tried changing the step size but to no avail. Isn't there any better method to find the Basin of attraction?

• Your question is not clear. If you try to reproduce results from the paper linked, then you should know that black hole horizon is given at r=2. So you should look any trajectory with initial condition r[0]>2. Mar 2 at 16:41
• @AlexTrounev, ah you are absolutely right. I made a blunder in the example. Let me edit that. My question is still the same though, how to effectively decide the initial values. Mar 3 at 1:45
• Could you remove wrong range {R, 1, 2, 0.01} from your post as well? Mar 3 at 2:08
• Ah yes sorry my mistake. Let me do that as well. Mar 3 at 2:50
• Thank you very much (+1). Mar 3 at 3:14

A somewhat dirty solution is to do the following:

tab = Quiet@
Table[If[Flatten[solution[R, 0.0]] === {}, Nothing,
solution[R, 0.0] // Last // First], {R, 2.72, 3.75, 0.01}]
• @codebpr Now I understand better what you are looking for, the basin of attraction is another task. :) Mar 2 at 7:27
• @E.Chan-López Please, pay attention that range {R, 1, 2, 0.01} is wrong in this problem since BH horizon is at R=2. Mar 3 at 2:11
• @codebpr Please ask about references to resource functions in a new post, as it would be useful for other users who wish to put resource functions to professional use. There is something additional I want to submit to update ClickPoincarePlot2D. Do you mind if I include you as a co-author? I'll take advantage of the fact that I'm going to include one of the systems that you considered in your posts, specifically the system where I added the solution using ClickPane. Mar 3 at 6:55
• It would be my pleasure to be a co-author, please include anything you like. I am also going to ask a new question regarding this now. Mar 3 at 7:06
• @codebpr I assume there is no LaTeX citation style requirement for resource functions. Surely the citation style is free and only requires adding the link in addition to the name of the authors. You can modify your question in this sense to avoid its closure. Mar 3 at 9:54

You may get the range directly from tab or the plot.

To get the range from "tab", you first flatten "tab" to get the points and then use MinMax:

tmp = Flatten[tab, 1];
MinMax /@ {tmp[[All, 1]], tmp[[All, 2]]}

{{2.23392, 3.91358}, {-26.4203, 167.055}}

To get the range from the plot: Say "pl" is your plot, then with "AbsoluteOptions":

AbsoluteOptions[p, PlotRange]

{PlotRange -> {{2.19893, 3.91358}, {-26.4203, 167.055}}}
• But that's when I already defined tab. My question is the other way around. How to find 'range' when I don't have a prior plot to look at? Basically, the basin of attraction. Mar 2 at 8:47
• I changed my answer and added the possibility to get the range from tab. Mar 2 at 9:54
• It works for this tab but sadly not for the modified one. Please see the edited question. Mar 2 at 13:04

Well, this problem is easily resolved by carefully choosing the approximate plot range using the ResourceFunction["ClickPoincarePlot2D"] developed by @E. Chan-López as shown below:

kr = 100; kθ = 25; rh = 2; rc = 3.2; θc = 0; κ =0.25;
f[r_] := 2  κ  (r - rh)

ham[r_, pr_, θ_, pθ_] := -Sqrt[1 - f[r]] pr + Sqrt[
pr^2 + pθ^2/r^2] + 1/2 kr (r - rc)^2 +
1/2 kθ (rh θ - rc θc)^2

H := Function[{S}, -S[[2]] Sqrt[1 - 0.5  (-2 + S[[1]])] + Sqrt[
S[[2]]^2 + S[[4]]^2/S[[1]]^2] + 50  (-3.2 + S[[1]])^2 +
25/2  (0 + 2  S[[3]])^2]

{a1, b1, c1,
d1} = {D[ham[r[t], pr[t], θ[t], pθ[t]], pr[t]],
D[ham[r[t], pr[t], θ[t], pθ[t]],
pθ[t]], -D[ham[r[t], pr[t], θ[t], pθ[t]],
r[t]], -D[
ham[r[t], pr[t], θ[t], pθ[t]], θ[
t]]}; hameqns = {r'[t] == a1, pr'[t] == c1, θ'[t] == b1,
pθ'[t] == d1};

cross = θ[t];
recover = {r[t], pr[t]};

ResourceFunction[
"ClickPoincarePlot2D"][hameqns, H, 75.0, t, 6000, cross, recover,
{PlotStyle -> AbsolutePointSize[1], PlotRange -> {{2, 4}, Automatic},
AspectRatio -> 1, PlotHighlighting -> None}]