Waveform Relaxation

The underlying idea behind wave form relaxation is to (as directly as possible implement the Picard-Lindelof iteration
	(Subscript[y, n+1])[t]=y[0]+(∫_(τ=0)^(τ=t))f((Subscript[y, n])[τ])dτ
Here is a direct implementation on our gravity problem with a constant initial condition.

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

Clear[zOld, t, z, zSol]
zOld[0][t_]=IC;
nIterMax=120
Do[
zSol=z/.NDSolve[{
z'[t]== Gravity[zOld[n][t]],
z[0]==IC
} ,
z,
{t, 0, TMax}
][[1]];
zOld[n+1]=zSol,
{n,0,nIterMax}]

Manipulate[
Show[
MmaPic,
ParametricPlot[
zOld[n][t][[{1,2}]],
{t,0, TMax},
PlotStyle->Red,
PlotPoints->50
]
],
{n, 1, 120,1}]

Ideas:
1) Do not iterate on the entire region at a time. 
2) Stop iterating on a chunk once it has converged.
3) Worry about how the integration in the "quadrature" is handled.
4) If we get fancy to adapt the size of the chunk we are using.

Here is an indication of a serious issue fr quite a few applications.

n1=nIterMax;
LogPlot[ 
Abs[ySol[t][[{1,2}]]-zOld[n1][t][[{1,2}]]],{t, 0, TMax},
AxesLabel->{"t","err"}
]

n1=nIterMax;
LogPlot[ 
Abs[ySol[t][[{1,2}]]-zOld[n1][t][[{1,2}]]],{t, 0, 0.001},
AxesLabel->{"t","err"}
]

Notation for Step size control

A) (Subscript[H, n]) is the size of the chunk we are iterating on at iteration n.
B) (Subscript[s, n]) is the start of the chunk we are iterating on at iteration n.
So what we have is that our iteration is 
	(Subscript[y, n+1])[t]=y[(Subscript[s, n])]+(∫_(τ=(Subscript[s, n]))^(τ=t))f((Subscript[y, n])[τ])dτ
unitll we advance.
Label the internal quadrature nodes somehow.