In this problem, we will simulate the motion of a planet orbiting the sun:

As shown, we place the sun is located at the origin, and the aphelion $($the point at which the planet is farthest away from the sun$)$ is located on the positive $x-$axis. At time $t,$ the planet is located at the coordinate $(x,y)$ and its distance from the sun is $r=\sqrt[\phantom{2}]{x^2+y^2}.$ We can model the motion of the planet as a set of two second$-$order ODEs:

$\frac{d^{2}x}{d{t}^2}=-GM(\frac{x}{r^3})$

$\frac{d^{2}y}{d{t}^2}=-GM(\frac{y}{r^3})$

where $G$ is the gravitational constant, $M$ is the mass of the sun, and the combination $GM$ is equal to $2\mathrm{.}94{10}^{-4}A{U}^3{d}^{-2}.($For this problem, we will use the time unit of days $($ d $),$ and the length unit of astronomical units $($ AU $),$ which is defined as the average distance from the Earth to the sun, about $149\mathrm{.}6{10}^{6}km.)$ We will choose the initial location of the planet to be the aphelion, namely,