The Henon-Heiles Hamiltonian, discussed in Lecture 34.  These commands can be used to produce all plots shown in class.

You'll need to run the "RootSearch.m" code, linked to on the class site.  It may throw an error (written for old version of Mathematica) but seems to work fine in practice.  When trying this code out for yourself, make sure to keep energy below 1/6 - above this value, the model is unstable.

A non-chaotic solution, with (presumably) quasiperiodic dynamics

x0=0.2;
y0=0.14;
p0=0.09;
q0=0.2;
tmax=805;
energy=0.5*(x0^2+y0^2+p0^2+q0^2)+x0^2*y0-y0^3/3
s=NDSolve[{x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]-2x[t]*y[t],q'[t]==-y[t]-x[t]^2+y[t]^2,p[0]==p0,q[0]==q0,x[0]==x0,y[0]==y0},{x,y,p,q},{t,0,tmax}];
ParametricPlot[Evaluate[{y[t],q[t]}/.s],{t,0,200},AxesLabel->{y,q}]
Plot[Evaluate[{y[t]}/.s],{t,0,tmax},AxesLabel->{t,y}]

A chaotic solution:

x0=0.2;
y0=0.44;
p0=0.09;
q0=0.2;
tmax=805;
energy=0.5*(x0^2+y0^2+p0^2+q0^2)+x0^2*y0-y0^3/3
s=NDSolve[{x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]-2x[t]*y[t],q'[t]==-y[t]-x[t]^2+y[t]^2,p[0]==p0,q[0]==q0,x[0]==x0,y[0]==y0},{x,y,p,q},{t,0,tmax}];
ParametricPlot[Evaluate[{y[t],q[t]}/.s],{t,0,200},AxesLabel->{y,q}]
Plot[Evaluate[{y[t]}/.s],{t,0,tmax},AxesLabel->{t,y}]

Here are “Poincare sections” of the dynamics, showing the transition from integrability to chaos upon increasing the energy of the system

x0=0.0;
y0=0.0;
p0=0.55*Sin[0.3];
q0=-0.55*Cos[0.3];
tmax=9000;
energy=0.5*(x0^2+y0^2+p0^2+q0^2)+x0^2*y0-y0^3/3
s=NDSolve[{x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]-2x[t]*y[t],q'[t]==-y[t]-x[t]^2+y[t]^2,p[0]==p0,q[0]==q0,x[0]==x0,y[0]==y0},{x,y,p,q},{t,0,tmax}];
xx[t_?NumericQ]:=(x[t]/.s)[[1]];
sol=RootSearch[xx[t]==0,{t,0,tmax}];
w=({y[t],q[t]}/.s)/.sol;
W1=ArrayReshape[w,{Length[w],2}];
x0=0.0;
y0=0.0;
p0=0.35*Sin[0.3];
q0=-0.35*Cos[0.3];
tmax=9000;
energy=0.5*(x0^2+y0^2+p0^2+q0^2)+x0^2*y0-y0^3/3
s=NDSolve[{x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]-2x[t]*y[t],q'[t]==-y[t]-x[t]^2+y[t]^2,p[0]==p0,q[0]==q0,x[0]==x0,y[0]==y0},{x,y,p,q},{t,0,tmax}];
xx[t_?NumericQ]:=(x[t]/.s)[[1]];
sol=RootSearch[xx[t]==0,{t,0,tmax}];
w=({y[t],q[t]}/.s)/.sol;
W2=ArrayReshape[w,{Length[w],2}];
x0=0.0;
y0=0.0;
p0=0.25*Sin[0.3];
q0=-0.25*Cos[0.3];
tmax=9000;
energy=0.5*(x0^2+y0^2+p0^2+q0^2)+x0^2*y0-y0^3/3
s=NDSolve[{x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]-2x[t]*y[t],q'[t]==-y[t]-x[t]^2+y[t]^2,p[0]==p0,q[0]==q0,x[0]==x0,y[0]==y0},{x,y,p,q},{t,0,tmax}];
xx[t_?NumericQ]:=(x[t]/.s)[[1]];
sol=RootSearch[xx[t]==0,{t,0,tmax}];
w=({y[t],q[t]}/.s)/.sol;
W3=ArrayReshape[w,{Length[w],2}];
ListPlot[{W1,W2,W3},PlotMarkers->Automatic,AxesLabel->{y,q},PlotStyle->{Red,Orange,Blue}]

Here are “Poincare sections” of the dynamics, showing the simultaneous existence of chaos and integrability at the same energy!

x0=0.0;
y0=0.0;
p0=0.5*Sin[0.3];
q0=-0.5*Cos[0.3];
tmax=9000;
energy=0.5*(x0^2+y0^2+p0^2+q0^2)+x0^2*y0-y0^3/3
s=NDSolve[{x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]-2x[t]*y[t],q'[t]==-y[t]-x[t]^2+y[t]^2,p[0]==p0,q[0]==q0,x[0]==x0,y[0]==y0},{x,y,p,q},{t,0,tmax}];
xx[t_?NumericQ]:=(x[t]/.s)[[1]];
sol=RootSearch[xx[t]==0,{t,0,tmax}];
w=({y[t],q[t]}/.s)/.sol;
W1=ArrayReshape[w,{Length[w],2}];
x0=0.0;
y0=0.0;
p0=0.5*Sin[0.09];
q0=-0.5*Cos[0.09];
tmax=9000;
energy=0.5*(x0^2+y0^2+p0^2+q0^2)+x0^2*y0-y0^3/3
s=NDSolve[{x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]-2x[t]*y[t],q'[t]==-y[t]-x[t]^2+y[t]^2,p[0]==p0,q[0]==q0,x[0]==x0,y[0]==y0},{x,y,p,q},{t,0,tmax}];
xx[t_?NumericQ]:=(x[t]/.s)[[1]];
sol=RootSearch[xx[t]==0,{t,0,tmax}];
w=({y[t],q[t]}/.s)/.sol;
W2=ArrayReshape[w,{Length[w],2}];
x0=0.0;
y0=0.0;
p0=0.5*Sin[0.5];
q0=-0.5*Cos[0.5];
tmax=9000;
energy=0.5*(x0^2+y0^2+p0^2+q0^2)+x0^2*y0-y0^3/3
s=NDSolve[{x'[t]==p[t],y'[t]==q[t],p'[t]==-x[t]-2x[t]*y[t],q'[t]==-y[t]-x[t]^2+y[t]^2,p[0]==p0,q[0]==q0,x[0]==x0,y[0]==y0},{x,y,p,q},{t,0,tmax}];
xx[t_?NumericQ]:=(x[t]/.s)[[1]];
sol=RootSearch[xx[t]==0,{t,0,tmax}];
w=({y[t],q[t]}/.s)/.sol;
W3=ArrayReshape[w,{Length[w],2}];
ListPlot[{W2,W1,W3},PlotMarkers->Automatic,AxesLabel->{y,q},PlotStyle->{Red,Orange,Blue}]