Euler and RK2

EulerStep[y_,{f_,h_}]:=y+h*f[y]
EulerSolveF[y0_,TMax_,{f_,n_}]:= Module[
{h=N[TMax/n],yy=y0},
Do[
yy= EulerStep[yy,{f,h}],
{n}];
yy
]
RK2Step[y_,{f_,h_}]:= Module[
{k1,k2},
k1=f[y];
k2=f[y+h k1];
y+(h/2)*(k1+k2)
]
RKSolveF[y0_,TMax_,{f_,n_}]:= Module[
{h=N[TMax/n],yy=y0},
Do[
yy= RK2Step[yy,{f,h}],
{n}];
yy
]

Gravitational Problem: p29 Euler vs RK FixedStep Size

f1[{y1_,y2_,y3_,y4_,y5_,y6_},μ_] :={y4,y5,y6,
2y5 + y1 - ((μ(y1+μ-1))/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)(y1+μ))/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2))),
-2y4 + y2 - ((μ y2)/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)y2)/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2))),
- ((μ y3)/((((y2^2)+(y3^2)+(((y1+μ-1))^2)))^(3/2))) -(((1-μ)y3)/((((y2^2)+(y3^2)+(((y1+μ))^2)))^(3/2)))
}
IC = {0.994,0,0,0,-2.001585106370825224,0};
TMax=17.06521656015796; 
μVal =1/81.45;
ySol[t_]=y[t]/.NDSolve[{
y'[t]==f1[y[t],μVal],
y[0]==IC
},
y[t],
{t, 0, TMax}
][[1]];

T = 0.01 TMax;
f[y_]=f1[y,μVal];
EulerData=Table[
n=(2^p);
yEuler=EulerSolveF[IC,T,{f, n}];
{T/n,Norm[ySol[T]-yEuler]
},
{p, 7, 15}]
RKData=Table[
n=(2^p);
yEuler=RKSolveF[IC,T,{f, n}];
{T/n,Norm[ySol[T]-yEuler]
},
{p, 7, 15}]
ListLogLogPlot[{EulerData,RKData},
Frame->True,
FrameLabel->{"h","error"},
PlotStyle->{Red, Blue},
GridLines->Automatic]

Stiffish Problem:  p60

f2[{t_,y1_,y2_}]:={1,-16 y1 + 12 y2 + 16 Cos[t] - 13 Sin[t],12 y1 - 9 y2 -11 Cos[t]+ 9 Sin[t]}
T =10.0;
IC={0.0,1.0,0.0};
f[y_]=f2[y];
ySol[t_]={t,Cos[t], Sin[t]};
EulerData=Table[
n=(1.5^p);
yEuler=EulerSolveF[IC,T,{f, n}];
{T/n,Norm[ySol[T]-yEuler]
},
{p, 4, 15}]
RKData=Table[
n=(1.5^p);
yEuler=RKSolveF[IC,T,{f, n}];
{T/n,Norm[ySol[T]-yEuler]
},
{p, 4, 15}]
ListLogLogPlot[{EulerData,RKData},
Frame->True,
PlotRange->All,
FrameLabel->{"h","error"},
PlotStyle->{Red, Blue},
GridLines->Automatic]

f2[{t_,y1_,y2_}]:={1,-16 y1 + 12 y2 + 16 Cos[t] - 13 Sin[t],12 y1 - 9 y2 -11 Cos[t]+ 9 Sin[t]}
T =0.01;
IC={0.0,1.0,0.0};
f[y_]=f2[y];
ySol[t_]={t,Cos[t], Sin[t]};
ParametricPlot3D[ ySol[t],{t, 0, 34},
BoxRatios->1]
EulerData=Table[
n=(1.5^p);
yEuler=EulerSolveF[IC,T,{f, n}];
{T/n,Norm[ySol[T]-yEuler]
},
{p, 4, 15}]
RKData=Table[
n=(1.5^p);
yEuler=RKSolveF[IC,T,{f, n}];
{T/n,Norm[ySol[T]-yEuler]
},
{p, 4, 15}]
ListLogLogPlot[{EulerData,RKData},
Frame->True,
PlotRange->All,
FrameLabel->{"h","error"},
PlotStyle->{Blue, Red},
GridLines->Automatic]