Preamble

Involves code snippets from:

Weisstein, Eric W. “Surface of Section.” From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/SurfaceofSection.html 
Enrique Zeleny “The Hénon-Heiles System”, http://demonstrations.wolfram.com/TheHenonHeilesSystem/

solveHenons[x0_,y0_,py_,e_,tf_]:=Reap@NDSolve[{{
(Subscript[y, 1])'[t]==(Subscript[y, 3])[t],
(Subscript[y, 2])'[t]==(Subscript[y, 4])[t],
(Subscript[y, 3])'[t]==-(Subscript[y, 1])[t] (1+2 (Subscript[y, 2])[t]),
(Subscript[y, 4])'[t]==-(((Subscript[y, 1])[t])^2)+(-1+(Subscript[y, 2])[t]) (Subscript[y, 2])[t]},
{(Subscript[y, 1])[0]==x0,(Subscript[y, 2])[0]==y0,(Subscript[y, 3])[0]==Sqrt[2/3 y0^3+2e-x0^2-y0^2-2x0^2y0-py^2],(Subscript[y, 4])[0]==py}},
{(Subscript[y, 1])[t],(Subscript[y, 3])[t],(Subscript[y, 2])[t],(Subscript[y, 4])[t]}, {t, tf},Method->{"EventLocator","Event"->(Subscript[y, 1])[t],"EventAction":> Sow[{(Subscript[y, 2])[t], (Subscript[y, 4])[t]}],"EventLocationMethod"->"LinearInterpolation",
"Method"->{"SymplecticPartitionedRungeKutta","DifferenceOrder"->4,"PositionVariables"->{(Subscript[y, 1])[t], (Subscript[y, 2])[t]}}},
StartingStepSize->0.25, MaxSteps->∞]

energyHenons[sdata_]:=(((Subscript[y, 3])[t])^2)/2+(((Subscript[y, 4])[t])^2)/2+(((Subscript[y, 1])[t])^2)/2+(((Subscript[y, 1])[t])^2)(Subscript[y, 2])[t]+(((Subscript[y, 2])[t])^2)/2-(((Subscript[y, 2])[t])^3)/3/.sdata[[1,1]]

showPoincare[sdata_]:=
If[sdata[[2]]=!={},ListPlot[sdata[[2,1]],Axes->False,Frame->True,AspectRatio->1,ImageSize->300,PlotStyle->{PointSize[.008]}]]

showOrbit[sdata_,tf_]:=ParametricPlot3D[Evaluate[{(Subscript[y, 1])[t],(Subscript[y, 2])[t], (Subscript[y, 4])[t]}/.First@sdata],{t,0,tf},
MaxRecursion -> 9,ViewPoint->{1,0,-6},Boxed->False,Axes->False,ImageSize->300,SphericalRegion-> True,ColorFunction->(ColorData["ThermometerColors"][#4]&)]

visualize[x0_,y0_,py_,e_,tf_]:=With[{sdata=solveHenons[x0,y0,py,e,tf]},
Row[{showOrbit[sdata,tf],showPoincare[sdata]}]
]

surfaceOfSectionEnergy[e_,tf_,npoints_]:=With[{x=.3},
Internal`DeactivateMessages[
ListPlot[Join@@Table[
With[{sdata=solveHenons[Random[Real,{-x,x}],Random[Real,{-x,x}],Random[Real,{-x,x}],e,tf]},
If[sdata[[2]]=!={},sdata[[2,1]],{}]],
{npoints}],
PlotStyle->{PointSize[.004]},
AspectRatio->1,AxesLabel->TraditionalForm/@{y[t],(Subscript[p, y])[t]},
ImageSize->500]
]
]

Henon-Heiles System

Manipulate[
visualize[x0,y0,py,e,tf],
{{tf,2000,Style["t", Italic]},4,50000, Appearance-> "Labeled"},
{{x0,0,Subscript[Style["x", Italic],"0"]},-.3,.3,Appearance-> "Labeled"},
{{y0,0,Subscript[Style["y", Italic],"0"]},-.3,.3,Appearance-> "Labeled"},
{{e,1/12,"E"},0,1/6,Appearance-> "Labeled"},
{{py,0,Subscript[Style["p", Italic],Style["y", Italic]]},-.3,.3,Appearance-> "Labeled"}]

visualize[0.02,0.02,0.1,0.02,5000]

visualize[0,0,0,1/12,20000]

visualize[0,0,0,0.088,2000]

Chaos emerging after long time.

visualize[0.02,0.02,0.19,1/8,20000]
visualize[0.02,0.02,0.19,1/8,30000]
visualize[0.02,0.02,0.19,1/8,100000]

Check that energy remains conserved despite numerical errors:

Plot[Evaluate[energyHenons[solveHenons[0.02,0.02,0.19,1/8,100000]]],{t,0,100000}]

Seems OK.

With[{x=0.0001},
Internal`DeactivateMessages[
ListPlot[Join@@Table[
With[{sdata=solveHenons[.02+Random[Real,{-x,x}],.02+Random[Real,{-x,x}],.19+Random[Real,{-x,x}],1/8,5000]},
If[sdata[[2]]=!={},sdata[[2,1]],{}]],
{100}],
PlotStyle->{PointSize[.004]},
AspectRatio->1,AxesLabel->TraditionalForm/@{y[t],(Subscript[p, y])[t]},
ImageSize->500]
]
]

visualize[0.02,0.02,0.19,1/6,1000]

visualize[0.02,0.02,0.19,1/6,100000]

surfaceOfSectionEnergy[1/12,5000,50]
surfaceOfSectionEnergy[1/8,5000,50]
surfaceOfSectionEnergy[1/6,5000,50]