With the following code, you can plot the Poincaré sections of the Hénon-Heiles system:

    With[{icv = {0, 0.36169437164930385`, 0.20100851639176504`, 0.029106357137938632`}}, 
    psection = Reap[NDSolve[{x'[t] == px[t], px'[t] == -(x[t] + 2 x[t]* y[t]), 
    y'[t] == py[t], py'[t] == -(y[t] + x[t]^2 - y[t]^2), 
    x[0] == icv[[1]], px[0] == icv[[2]], y[0] == icv[[3]], 
    py[0] == icv[[4]]}, {x, px, y, py}, {t, 0, 1000}, 
    MaxSteps -> ∞, 
    Method -> {"EventLocator", "Event" -> x[t], 
    "EventAction" :> Sow[{y[t], py[t]}]}]][[2]];]

We load the MaTeX package for labels with LaTeX:

    Needs["MaTeX`"];

Finally, we plot de Poincaré section:

    ListPlot[section, PlotRange -> All, AspectRatio -> 1, PlotStyle -> Black, 
    Frame -> True, FrameStyle -> Black, Axes -> False, LabelStyle -> Directive[Black, Small], 
    FrameLabel -> {{MaTeX["p_{y}", Magnification -> 1.5], None}, 
    {MaTeX["y", Magnification -> 1.5], 
    MaTeX["\\text{Hénon-Heiles system}", Magnification -> 1.3]}}, 
    RotateLabel -> False, Epilog -> Inset[MaTeX["E=0.08333", Magnification -> 1], 
    {0.21, 0.15}, Automatic, 1], ImageSize -> Medium]

The Poincaré section:

[![enter image description here][1]][1]

You can use a For loop to compute with more initial conditions subject to the energy constraint of the Hénon-Heiles system. For more details, see my answer [Poincaré Sections for spring pendulum][2].

With many initial conditions:

    With[{icv = {{0, 0.4082401254164024`, 0, 0}, 
    {0, 0.24008038662195985`, -0.15316368171208566`, 0.28838672839135326`}, 
    {0, 0.24352018790129148`, 0.09848498536170303`, 0.3135210500944904`}, 
    {0, 0.20788819727584018`, 0.18687849877360435`, 0.3047456328865395`}, 
    {0, 0.11124773284074371`, 0.24234454959281115`, 0.32410152798705066`}, 
    {0, 0.09612790055948518`, 0.2562151235263743`, 0.32091473672760257`}, 
    {0, 0.14504873906640348`, 0.2433184249171519`, 0.3098719090788399`}, 
    {0, 0.21466692008332666`, 0.21341797098009158`, 0.2855018063099425`}, 
    {0, 0.2242200636897092`, 0.21057030870037158`, 0.27976766448196655`}, 
    {0, 0.3683386516212832`, 0.18778901046261137`, 0.011696534024095057`}, 
    {0, 0.3451153373835715`, 0.23762310035771178`, 0.005962392196119196`}, 
    {0, 0.33504538978074816`, 0.256132905175892`, 0.0016617858251372995`}, 
    {0, 0.37097564355830726`, -0.03575247849541091`, 0.1665183633794433`}, 
    {0, 0.37507484064199426`, -0.06280527015275122`, 0.14788240243852177`}, 
    {0, 0.3774580162329546`, -0.08416273725065143`, 0.12924644149760023`}, 
    {0, 0.3779249463635693`, -0.09982487978911159`, 0.11491108692766056`}, 
    {0, 0.3818158505457779`, -0.10979169776813169`, 0.08910744870176919`}, 
    {0, 0.3850735819245613`, -0.11406319118771174`, 0.06617088138986574`}, 
    {0, 0.3571135001050445`, -0.02151416709681081`, 0.1966226079763166`}, 
    {0, 0.3537329768656702`, -0.0015805311387706023`, 0.20379028526128642`}}}, 
    psection = Reap[
    Table[NDSolve[{x'[t] == px[t], px'[t] == -(x[t] + 2 x[t]* y[t]), 
    y'[t] == py[t], py'[t] == -(y[t] + x[t]^2 - y[t]^2), 
    x[0] == Part[icv[[i]], 1], px[0] == Part[icv[[i]], 2], 
    y[0] == Part[icv[[i]], 3], py[0] == Part[icv[[i]], 4]}, {x, 
    px, y, py}, {t, 0, 5000}, MaxSteps -> \[Infinity], 
    Method -> {"EventLocator", "Event" -> x[t], 
    "EventAction" :> Sow[{y[t], py[t]}]}], {i, 1, Length[icv]}]][[2]];]

[![enter image description here][3]][3]

Trajectories:

    With[{icv = {{0, 0.4082401254164024`, 0, 0}, 
    {0, 0.24008038662195985`, -0.15316368171208566`, 0.28838672839135326`}, 
    {0, 0.24352018790129148`, 0.09848498536170303`, 0.3135210500944904`}, 
    {0, 0.20788819727584018`, 0.18687849877360435`, 0.3047456328865395`}, 
    {0, 0.11124773284074371`, 0.24234454959281115`, 0.32410152798705066`}, 
    {0, 0.09612790055948518`, 0.2562151235263743`, 0.32091473672760257`}, 
    {0, 0.14504873906640348`, 0.2433184249171519`, 0.3098719090788399`}, 
    {0, 0.21466692008332666`, 0.21341797098009158`, 0.2855018063099425`}, 
    {0, 0.2242200636897092`, 0.21057030870037158`, 0.27976766448196655`}, 
    {0, 0.3683386516212832`, 0.18778901046261137`, 0.011696534024095057`}, 
    {0, 0.3451153373835715`, 0.23762310035771178`, 0.005962392196119196`}, 
    {0, 0.33504538978074816`, 0.256132905175892`, 0.0016617858251372995`}, 
    {0, 0.37097564355830726`, -0.03575247849541091`, 0.1665183633794433`}, 
    {0, 0.37507484064199426`, -0.06280527015275122`, 0.14788240243852177`}, 
    {0, 0.3774580162329546`, -0.08416273725065143`, 0.12924644149760023`}, 
    {0, 0.3779249463635693`, -0.09982487978911159`, 0.11491108692766056`}, 
    {0, 0.3818158505457779`, -0.10979169776813169`, 0.08910744870176919`}, 
    {0, 0.3850735819245613`, -0.11406319118771174`, 0.06617088138986574`}, 
    {0, 0.3571135001050445`, -0.02151416709681081`, 0.1966226079763166`}, 
    {0, 0.3537329768656702`, -0.0015805311387706023`, 0.20379028526128642`}}}, 
    trajectories = Reap[
    Table[NDSolve[{x'[t] == px[t], px'[t] == -(x[t] + 2 x[t]* y[t]), 
    y'[t] == py[t], py'[t] == -(y[t] + x[t]^2 - y[t]^2), 
    x[0] == Part[icv[[i]], 1], px[0] == Part[icv[[i]], 2], 
    y[0] == Part[icv[[i]], 3], py[0] == Part[icv[[i]], 4]}, {x, 
    px, y, py}, {t, 0, 5000}, MaxSteps -> \[Infinity], 
    Method -> {"EventLocator", "Event" -> x[t], 
    "EventAction" :> Sow[{y[t], py[t]}]}], {i, 1, Length[icv]}]][[1]];]
    (*A trajectory*)
    tr = ParametricPlot[Evaluate[{x[t], y[t]} /. trajectory[[2]]], {t, 0, 500}, PlotPoints -> 500, PlotStyle -> {Blue,Thickness[0.002]}, AxesStyle -> Black, LabelStyle -> Directive[Black, Small]];

A trajectory:

[![enter image description here][4]][4]

Chaotic regime:

[![enter image description here][5]][5]

Interactive way:

[![enter image description here][6]][6]


  [1]: https://i.sstatic.net/YbCmd.png
  [2]: https://mathematica.stackexchange.com/questions/61637/poincare-section-of-an-hamiltonian/241621#241621
  [3]: https://i.sstatic.net/SYr4D.png
  [4]: https://i.sstatic.net/XbJdA.jpg
  [5]: https://i.sstatic.net/DhxAM.png
  [6]: https://i.sstatic.net/nRPLj.gif