tmax=35  (* time of run *)
Z = 2  (* Let's do Helium *)

(* The Hamiltonian *)
H = (x2[t]^2 p1[t]^2 + x1[t]^2 p2[t]^2)/8 - Z (x1[t]^2+x2[t]^2) +
 	(x1[t] x2[t])^2 (1+ 1/(x1[t]^2+x2[t]^2))

(* Hamilton's Equations *)
HE = {x1'[t] == D[H, p1[t]], x2'[t] == D[H, p2[t]], 
	p1'[t] == -D[H, x1[t]], p2'[t] == -D[H, x2[t]]}

(* Initial Conditions *)
IC = {x1[0] == 2.0, x2[0]==-1.0, p1[0]==2.0, H==0 /.t->0}

(* The Solution *)
Solution = NDSolve[Join[HE,IC], {x1[t], x2[t], p1[t], p2[t]},
	{t, 0, tmax}, MaxSteps->10000]

(* Four Graphs *)
(* x1^2 versus x2^2 *)
ParametricPlot[Evaluate[{x1[t]^2, x2[t]^2} /. Solution], {t,0,tmax}]

(* The trajectory projected into the space x1, p1, x2 *)
Poincare = ParametricPlot3D[Evaluate[{x1[t], p1[t], x2[t]} /. Solution],
	 {t,0,tmax}]

(* A Poincare Section of the Trajectory *)
Show[Poincare, PlotRange->{-0.01,0.01}]

(* Check the Hamiltonian to see how well the integration program does *)
Plot[Evaluate[H /. Solution], {t,0,tmax}, PlotPoints->5,PlotDivision->10]