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]]}

(* Solve Hamilton's equations starting from x1=x2=x and with p1=-p2 *)
SolveHelium[x_, tmax_] := Block[{IC},
  IC = Solve[{x1[0] ==  x2[0] == x , p1[0]== -p2[0], H==0 /.t->0}][[2]] /. 
                                                               Rule->Equal;
  NDSolve[Join[HE,IC], {x1[t], x2[t], p1[t], p2[t]},
	{t, 0, tmax}]
]


Off[ParametricPlot::ppcom] (* switch off annoying error message *)

(* Plot trajectory in x1^2 and x2^2 starting from x1=x2=x and p1=-p2 *)
PlotHelium[x_, tmax_] := Block[{Solution},
  Solution = SolveHelium[x, tmax];
  ParametricPlot[Evaluate[{1.1 x^2 {t,t}/tmax,
                           {x1[t]^2, x2[t]^2} /. Solution}],
                           {t,0,tmax}, AspectRatio->1,
                           PlotRange->All,
                           PlotLabel->{"x= ", x, " t= ", tmax}]
]

PlotHelium[l_List] := PlotHelium[l[[1]], l[[2]]]

(* Given a good guess for initial and final period of the orbit plus a   *
 * time longer than the periodic orbit and an improvement factor returns *
 * a better guess provided the factor is small enough                    *)
Clean[{xguess_, tguess_}, tmax_, factor_] := Block[{ht, rt, K, t},
  Solution = SolveHelium[xguess, tmax];
  K[t_] := Evaluate[x1[t]^2 - x2[t]^2 /. Solution][[1]];
  rt = FindRoot[K[t]==0, {t, tguess}][[1]];
  {xguess - factor*(xguess - ((Abs[x1[t]] /. Solution) /. rt)[[1]]), t /. rt}
]

(* Iteratively applies Clean to find a better estimate for the starting *
 * position and period of the periodic orbit.  Plots a graph of the     *
 * new guesses to check for convergence                                 *)
CleanRoot[xguess_, tguess_, tmax_, factor_, n_] := Block[{list},
  list = NestList[Clean[#, tmax, factor]&, {xguess, tguess}, n];
  ListPlot[Transpose[list][[1]], Prolog->AbsolutePointSize[5]];
  Print[N[list,10]];
  N[list[[-1]], 10]
]

(* Given an initial and starting position, x,  integrates the action up *
 * to time tmax                                                         *)
Action[x_, tmax_] := Block[{po, t},
  po = SolveHelium[x, tmax][[1]];
  NIntegrate[x1[t]*x1[t]*x2[t]*x2[t]/Pi /. po, {t,0,tmax}]
]

Action[l_List] := Action[l[[1]], l[[2]]]