SV[g_,{y0_, v0_},{h_,n_}]:= Module[
{yData=ConstantArray[0,{n}],i},
yData[[1]]=y0;
yData[[2]]=y0 + h v0+ (1/2)(h^2) g[y0];
Do[
yData[[i+1]]=2yData[[i]]-yData[[i-1]] + (h^2)g[yData[[i]]],
{i, 2, n-1}];
yData
]

Euler[f_,y0_,{h_,n_}]:= Module[
{yData=ConstantArray[0,{n}],i},
yData[[1]]=y0;
Do[
yData[[i+1]]=yData[[i]] + h f[yData[[i]]],
{i, 1, n-1}];
yData[[All,{1,2}]]
]

RK4[f_,y0_,{h_,n_}]:= Module[
{yData=ConstantArray[0,{n}],f1,f2,f3,f4,i},
yData[[1]]=y0;
Do[
f1=f[yData[[i]]];
f2=f[yData[[i]]+(1/2)h f1];
f3=f[yData[[i]]+(1/2) h f2];
f4=f[yData[[i]]+h  f3];
yData[[i+1]]=yData[[i]] + (h/6) (f1+2f2+2f3+f4),
{i, 1, n-1}];
yData[[All,{1,2}]]
]

f[{y1_,y2_,y3_,y4_}]:= Module[
{d=((((y1^2)+(y2^2)))^(-3/2))},
{y3,y4,-y1 d, -y2 d}
]
g[{y1_,y2_}]:= Module[
{d=((((y1^2)+(y2^2)))^(-3/2))},
d{-y1, -y2}
]
y1st[e_]:= {1-e,0,0,(Sqrt[((1+e)/(1-e))])}
y0[e_]:= {1-e,0}
v0[e_]:={0,(Sqrt[((1+e)/(1-e))])}
e=0.1;
h=0.2;
TMax= 20*2π;
n=Floor[(TMax/h)];
SVData=SV[ g, {y0[e],v0[e]},{h, n}];
EulerData=Euler[ f, y1st[e],{h, n}];
RK4Data=RK4[ f, y1st[e],{h, n}];
ListPlot[{SVData,RK4Data},
Joined->True,
PlotStyle->{Red, Blue, Green}]

Now clearly the RK4 and the SV are doing a much better job than the Euler!

ListPlot[ {SVData[[All,1]],RK4Data[[All,1]]},
PlotStyle->{Red, Blue,Green}
]
ListPlot[ {SVData[[All,2]],RK4Data[[All,2]]},
PlotStyle->{Red, Blue, Green}
]

The system is supposed to conserve the energy
	(1/2)||v((||)^2)-(( 1)/(Norm[y]))
and we can see some interesting things if we make functions and look plot the result.

ySV[t_]=Map[ ListInterpolation[#,{{0, TMax}}][t]&, Transpose[SVData]];
yRK4[t_]=Map[ ListInterpolation[#,{{0, TMax}}][t]&, Transpose[RK4Data]];
Plot[{
(1/2)((Norm[v0[e]])^2)-(1/(Norm[y0[e]])),
(1/2)((Norm[ySV'[t]])^2)-(1/(Norm[ySV[t]])),
(1/2)((Norm[yRK4'[t]])^2)-(1/(Norm[yRK4[t]]))
},
{t, 0, TMax},
PlotStyle->{Red, Blue,Green},
PlotRange->All,
Frame->True,
Axes->False]

Apparently the initial conditions are slightly different and in addition the SV curve although it bounces around is not drifting.