我试图使用matlab绘制水星近日点的进动。为此，我正在关注Nicholas J. Giordano 和 Hisao Nakanishi 第二版的《计算物理学》一书。在那本书中，作者使用了广义相对论预测的以下力定律：

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

其中，它们用光速、太阳质量、轨道离心率来表示。$\alpha$

为了绘制轨道，我在 matlab 中使用了该方程和 ode45、ode23 函数。我的代码如下：

% 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])

ode45 或 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

使用 ode45 时，我得到以下情节：

如果我使用 ode23 我得到下图：

这些不是我所期待的情节。我不明白我的代码有什么问题。是我的代码还是 ODE 求解器？