The techniques described in this section are not for differential equations
containing nonsmooth functions. For example, you might have a differential
equation whose right-hand side involves a function that is evaluated by table look-up
and interpolation. If so, go back to Runge-Kutta with adaptive stepsize choice:
That method does an excellent job of feeling its way through rocky or discontinuous
terrain. It is also an excellent choice for quick-and-dirty, low-accuracy solution
of a set of equations. A second warning is that the techniques in this section are
not particularly good for differential equations that have singular points inside the
interval of integration. A regular solution must tiptoe very carefully across such
points. Runge-Kutta with adaptive stepsize can sometimes effect this; more generally,
there are special techniques available for such problems, beyond our scope here.
Apart from those two caveats, we believe that the Bulirsch-Stoer method,
discussed in this section, is the best known way to obtain high-accuracy solutions
to ordinary differential equations with minimal computational effort. (A possible
exception, infrequently encountered in practice, is discussed in §16.7.)
Numerical Recipes in C p718.

Modified Mid Point aka Gragg

ModifiedMidPointStep[yIn_,f_,{n_,H_}]:=Module[
{Data=ConstantArray[0,n+1],h=H/n,m},
Data[[1]]=yIn;
Data[[2]]=yIn+ h f[yIn];
Do[
Data[[m+1]]=Data[[m-1]]+ 2 h f[Data[[m]]],
{m,2,n}];
(1/2)(Data[[n+1]]+Data[[n]]+h f[Data[[n+1]]])
]

ModifiedMidptExtrapolation[yIn_,f_,{nSequence_List,H_}]:=Module[
{ExtrapolationData,n,i,hSq,vars},
ExtrapolationData = Table[ 
n=nSequence[[i]];
{((((H/n)))^2),ModifiedMidPointStep[yIn,f,{n,H}]},
{i, Length[nSequence]}
];
vars=Table[(hSq^(i-1)),{i,Length[nSequence]}];
Table[
Fit[
Transpose[{ExtrapolationData[[All,1]],ExtrapolationData[[All,2,i]]}],
vars,hSq],
{i, Length[yIn]}]/.hSq->0
]

p29 Test problem

f[{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.06521656015796; 
μVal =1/81.45;
ySol[t_]=y[t]/.NDSolve[{
y'[t]==f[y[t]],
y[0]==IC
},
y[t],
{t, 0, TMax},
MaxSteps->(10^5),
PrecisionGoal->12,
AccuracyGoal->12
][[1]];
ParametricPlot[ ySol[t][[{1,2}]],{t,0,TMax},
AxesLabel->{"y1","y2"}]

LTE

y=IC;
ConvData=Table[
{nInt,H}={(2^p),0.1};
{(H/nInt),Norm[ySol[H]-ModifiedMidPointStep[y,f,{nInt,H}]]},
{p,5,12}];
ListLogLogPlot[ConvData,
Frame->True,
FrameLabel->{"h","abs err"},
PlotLabel->"Modified Mid Point LTE"
]
Fit[Log[ConvData],{1,x},x]

Looks to me as though the LTE~(h^2) which we could verify analytically.  What is special about this scheme is that the LTE is a function of (h^2) not h.  The way you show this is through a time reversal argument.  So
	LTE = (Subscript[C, 1])(h^2)+(Subscript[C, 2])(h^4)+(Subscript[C, 3])(h^6)+…+(Subscript[C, i])(h^(2i))+…
and if our solution is smooth enough we can use solutions with different hs to solve and eleminate the Cs.  Doing this gives a particular extrapolation scheme called Bulrisch-Stoer which is probably the method of choice for extremely high accuracy

The techniques described in this section are not for differential equations
containing nonsmooth functions. For example, you might have a differential
equation whose right-hand side involves a function that is evaluated by table look-up
and interpolation. If so, go back to Runge-Kutta with adaptive stepsize choice:
That method does an excellent job of feeling its way through rocky or discontinuous
terrain. It is also an excellent choice for quick-and-dirty, low-accuracy solution
of a set of equations. A second warning is that the techniques in this section are
not particularly good for differential equations that have singular points inside the
interval of integration. A regular solution must tiptoe very carefully across such
points. Runge-Kutta with adaptive stepsize can sometimes effect this; more generally,
there are special techniques available for such problems, beyond our scope here.
Apart from those two caveats, we believe that the Bulirsch-Stoer method,
discussed in this section, is the best known way to obtain high-accuracy solutions
to ordinary differential equations with minimal computational effort. (A possible
exception, infrequently encountered in practice, is discussed in §16.7.)
Numerical Recipes in C p718.

Extrapolation

TStart=1.2;
y=ySol[TStart];
H=0.5;
nSequence={4,8,10,12};
ModifiedMidPointStep[y,f,{Max[nSequence],H}]
ModifiedMidptExtrapolation[y,f,{nSequence,H}]
ySol[TStart+H]