Code snippets for ODE Driver/Stepper

ODE23Step[y_,f_,h_]:= Module[
{k1,k2,k3,k4,yNew,zNew},
k1=f[y];
k2=f[y+ 1/2 h k1];
k3=f[y+3/4 h k2];
yNew=y + h ( 2/9 k1 + 1/3 k2 + 4/9 k3);
k4=f[yNew];
zNew=y +h*(7/24 k1 + 1/4 k2 + 1/3 k3 + 1/8 k4);
{yNew, Norm[zNew-yNew]}
]  

DriverFixed[y0_,f_,n_, TMax_,Stepper_]:=Module[
{Data,y=y0,err,h=TMax/n,t=0.0},
Data = Table[
{y, err} = Stepper[y, f, h];
t = t + h;
{t, y},
{n}
];
(* fixes it so that it comes out in the form {ts, ys} *)
Transpose[Data]
]

DriverVariableStep1[y0_,h0_,f_,n_, ErrorOrder_,LocalTol_,Stepper_]:=Module[
{y=y0,t=0.0,h=h0,err,Data, SafetyFactor=0.9, yNew},
  Data = Table[
{yNew, err} = Stepper[y, f, h];
If[err<LocalTol, {t,y}={t+h,yNew}];
h*=SafetyFactor*((((LocalTol/err)))^(1/ErrorOrder));
{t, y},
{n}];
(* fixes it so that it comes out in the form {ts, ys} *)
Transpose[Data]
]

Testing on our gravitaional problem

Testing this code on the gravitaional problem

f1[{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};

n=20;
TMax=0.5;
h0=TMax/n;
LocalTol=(10^(-4));
ErrorOrder = 3;
{tFs, yFs}=DriverFixed[IC,f1,n,TMax,ODE23Step];
{tVs, yVs}=DriverVariableStep1[IC,h0,f1,n,4,LocalTol,ODE23Step];
{ListPlot[{yFs[[All,{1,2}]],yVs[[All,{1,2}]]},
PlotStyle->{Red, Green},
PlotRange->All,
AxesLabel->{"y1","y2"}],
ListPlot[{tFs,tVs},
PlotStyle->{Red, Green},
PlotRange->All,
AxesLabel->{"i","t"}]}