WR-SIMP (Fixed Size Window)

The gravitational problem

IC={0.994,0.0,0.0,0.0,-2.0015851063790825224,0.0};
tmin=0.0;
tmaxold=17.06521656015796;
f[{y1_,y2_,y3_,y4_,y5_,y6_}]:=Module[{μ=1/81.45},
{
y4,
y5,
0,
2y5+y1-((μ(y1+μ-1))/(((y2^2+(((y1+μ-1))^2)))^(3/2)))-(((1-μ)(y1+μ))/(((y2^2+(((y1+μ))^2)))^(3/2))),
-2y4+y2-((μ*y2)/(((y2^2+(((y1+μ-1))^2)))^(3/2)))-(((1-μ)y2)/(((y2^2+(((y1+μ))^2)))^(3/2))),
0
}
];

The matrices: A and B

matrixA[r_]:=Module[{A,i,j},
A=ConstantArray[0,{r,r}];
For[i=2,i≤r-1,i++,
A[[i,1]]=1.0;
For[j=2,j≤i,j++,
A[[i+1,j]]=2.0;
];
A[[i,i]]=1.0;
];
A[[r,1]]=A[[r,r]]=1.0;
Return[A]
]

matrixB[r_,p_]:=Module[{A,i,j},
A=ConstantArray[0,{r,r-1}];
For[i=2,i≤r,i++,
For[j=2,j≤i,j++,
A[[i,j-1]]=1.0*p;
];
];
Return[A]
]

WR using Simpson' s rule

WRSIMP[r_,y_,e_]:=Module[{A,B,Yold,Ynew,F,Y0,h,a,g,k},
A=matrixA[r];
B=matrixB[r,4];
Y0= ConstantArray[y,{r}];
Ynew=Y0;
h=(1.0*e)/(2*r-2);
For[k=1,k≤100,k++,
F=Table[f[Ynew[[i]]],{i,1,r}];
g=Table[f[((Ynew[[i]]+Ynew[[i+1]])/2.0)+((h*(f[Ynew[[i]]]-f[Ynew[[i+1]]]))/4.0)],{i,1,r-1}];
Yold=Ynew;
Ynew=Y0+(h/3.0)*(A.F+B.g);
If[Norm[Yold-Ynew]≤(10^(-5)),Goto[endfor]];
];
Label[endfor];
Return[Ynew]
]

Mathematica NDSolve

MmaSol=y/.NDSolve[ {
y'[t]==f[y[t]],
y[0]==IC
},
y,
{t, 0, tmaxold}][[1]];
MmaPic=ParametricPlot[MmaSol[t][[{1,2}]],{t,0,tmaxold}];

Run on the test problem

Example1=Manipulate[
yin=MmaSol[tmin];
eval=Max[Abs[Eigenvalues[D[f[{y1,y2,y3,y4,y5,y6}],{{y1,y2,y3,y4,y5,y6}}]/.y1->yin[[1]]/.y2->yin[[2]]]]];
myPlot=ListLinePlot[WRSIMP[nodes,yin,H][[All,1;;2]],PlotStyle->{Red,PointSize[Large]},Mesh->All];
Grid[{{Show[MmaPic,myPlot]," "},{"","Max Eigenvalue:",eval}}],
{nodes,10,100,2},
{tmin,0,tmaxold},
{H,0.1,tmaxold}
];

Example1