Can't plot correctly precession of perihelion of Mercury in MATLAB using ode45 or ode23

Question

I was trying to plot precession of perihelion of Mercury using matlab. For this I am following a book Computational Physics by Nicholas J. Giordano and Hisao Nakanishi 2nd Edition. In that book authors used following force law predicted by general relativity:

$$ F \approx \frac{G M_m M_s}{r^2}\left(1+\frac{\alpha}{r^2}\right) $$

Where $\alpha$, they state to be expressed by speed of light, mass of sun, eccentricity of the orbit.

To plot the orbit I have used that equation and ode45, ode23 functions in matlab. My code is as follows:

% Initial Conditions
% Initial position in AU
xm0 = [0.47 0];
% Initial velocity in AU/Year
vm0 = [0 8.2];

inivec = [xm0;vm0]

% tspan for 1000 years
tspan = [0, 1000];

[t, y] = ode45(@sunmer, tspan, inivec);

% x positions of Mercury
xmer = y(:,[1:2]);
% y positions of Mercury
vmer = y(:,[3:4]);
plot(xmer(:,1),xmer(:,2))
hold on

% Position of Sun in plot
plot(0,0,'o')
hold off
axis([-0.5, 0.5, -0.5, 0.5])

function for ode45 or ode23

function xpr = sunmer(t,x)
alpha = 0.0008;
xmer = x([1:2]);
vmer = x([3:4]);
r = norm(xmer);
rcubed = r^3;
rsqr = r^2;

% amer = -((G*Ms*x)/r^3)*(1+(alpha/r^2))
% G*Ms = 4*pi^2 in AU^3/Yr^2
amer = -(((4*(pi)^2))/rcubed)*xmer*(1+(alpha/rsqr));
xpr = [vmer;amer];
end

I am getting following plot when using ode45:

If I use ode23 I get following figure:

These are not the plots I was expecting. I can't understand what is wrong in my code. Is it my code or the ODE solvers?