Conservation

Lots of ODE systems 
	y'[t]=f[y[t]]
folks are interested in have conserved properties.  In particular anything Newtonish
	f(y)= m y'' 
without damping has a conserved energy. Dot the equation with y' and integrate
	f(y).y'= m y'.y''
integrate 
	(d/dt)[F(y)]= (d/dt)[(1/2)m (((y'))^2)] 
to find
	(1/2)m (v^2)-F[y]=constant
for an exact solution.  A couple of good questions are: 
1) How well do our existing schemes do at preserving this energy?
2) Can I design special schemes to preserve this energy and other conserved quantities?

Back to our gravity test problem.

Gravity[{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{μ=1/81.45,d1,d2 },
d1 =((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2)); 
d2=((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2));
{y4,y5,y6,
2y5 + y1 - ((μ(y1+μ-1))/d1) -(((1-μ)(y1+μ))/d2),
-2y4 + y2 - ((μ y2)/d1) -(((1-μ)y2)/d2),
- ((μ y3)/d1) -(((1-μ)y3)/d2)
}
]
IC = {0.994,0,0,0,-2.001585106370825224,0};
TMax=3*17.5;
ySol=y/.NDSolve[ {
y'[t]==Gravity[y[t]],
y[0]==IC
},
y,
{t, 0, TMax}][[1]];
MmaPic=ParametricPlot[ySol[t][[{1,2}]],{t,0,TMax}]

In this case the conserved quantity is the formula below.

ConsGravity[{y1_,y2_,y3_,y4_,y5_,y6_}] :=Module[
{μ=1/81.45,d1,d2 },
d1 =((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(1/2)); 
d2=((((y2^2)+(y3^2)+(((y1+μ))^2)))^(1/2));
(1/2)((y4^2)+(y5^2)+(y6^2) -(y1^2)-(y2^2))- (μ/d1)-(((1-μ))/d2)
]
Plot[1.42821+ConsGravity[ySol[t]],{t, 0,TMax},
PlotRange->All]

This looks pretty well conserved. We should see how well our solvers do at conserving this "energy".