<s>[INST] Using the Runge Kutta's Method to solve a 2nd derivative question Given that the motion of two bodies subject to a gravitational force of $\frac{d^2x}{dt^2}$ = $\frac{-x}{(x^2+y^2)^{3/2}}$, $x(0) = 0$, $\frac{dx}{dt}$$(0) = -0.5$  and $\frac{d^2y}{dt^2}$ = $\frac{-y}{(x^2+y^2)^{3/2}}$, $y(0) = 1$, $\frac{dy}{dt}(0) = -0.5$ 
It is also given that at those points given, the motion is periodic. 
Find the largest and smallest $\sqrt{x^2+y^2}$(distance between the two bodies).
First thoughts: I am only experienced working with just first derivative so I'm not really sure if I am supposed to use the Runge Kutta method two times to find the original. I will also be computing later via matlab and not by hand as the computations can get extremely difficult. However, I am stuck on just where to begin. Or is creating a systems of solutions the proper way to go about this question?
EDIT: searching around the internet, I read that it might be possible to solve this using a system of equations? So I created two matrices:
//sorry I'm unaware of how to create matrices:
Matrix 1: X' = [$x_2$ $\frac{-x_1}{(x_1^2+y^2)^{3/2}}$] but this should be 1 columned and two rows.
Matrix 2: X(0) = $[0, -0.5]^T$
 [/INST] There is an alternative method to my previous answer when the acceleration function is defined as a 2nd order diff. equation $$\ddot{\mathbf{Y}} = f(t,\mathbf{Y},\dot{\mathbf{Y}})$$
with the state vector $\mathbf{Y} = (x,y)$ in your case, holding the positions in the two coordinates, and $\mathbf{V} = \dot{\mathbf{Y}}= (\dot{x},\dot{y})$ holding the velocity vectors
$$ \begin{pmatrix}\ddot{x}\\ \ddot{y}\end{pmatrix} = f(t,\begin{pmatrix}x\\y\end{pmatrix},\begin{pmatrix}\dot{x}\\ \dot{y}\end{pmatrix}) = \begin{pmatrix} -\frac{x}{\sqrt{(x^2+y^2)^3}} \\ -\frac{y}{\sqrt{(x^2+y^2)^3}}\end{pmatrix} $$
Then you follow the following steps.
$$\begin{aligned}
  \mathbf{C}_0 &= \mathbf{V} & \mathbf{K}_0 & = f(t, \mathbf{Y}, \mathbf{C}_0) \\
  \mathbf{C}_1 &= \mathbf{V}+ \frac{h}{2} \mathbf{K}_0 & \mathbf{K}_1 & = f(t+\frac{h}{2}, \mathbf{Y} + \frac{h}{2} \mathbf{C}_0, \mathbf{C}_1)\\
  \mathbf{C}_2 &= \mathbf{V}+ \frac{h}{2} \mathbf{K}_1 & \mathbf{K}_2 & = f(t+\frac{h}{2}, \mathbf{Y} + \frac{h}{2} \mathbf{C}_1, \mathbf{C}_2) \\
  \mathbf{C}_3 &= \mathbf{V}+ h \mathbf{K}_2 & \mathbf{K}_3 & = f(t+h, \mathbf{Y} + h \mathbf{C}_2, \mathbf{C}_3)
\end{aligned} $$
and
$$\begin{aligned}
\Delta \mathbf{Y} & = \frac{h}{6} \left( \mathbf{C}_0 + 2 \mathbf{C}_1 + 2 \mathbf{C}_2 + \mathbf{C}_3 \right) \\ 
\Delta \mathbf{V} & = \frac{h}{6} \left( \mathbf{K}_0 + 2 \mathbf{K}_1 + 2 \mathbf{K}_2 + \mathbf{K}_3 \right)
\end{aligned}$$
Essentially your are running two Rk4 schemes in parallel, one for the positions and one for the velocities and you keep all cross term influences appropriately.
NOTE: The position step can be simplified to $\Delta \mathbf{Y}  = h \left(\mathbf{V}+\frac{h}{6} \left( \mathbf{K}_0 +  \mathbf{K}_1 +  \mathbf{K}_2  \right)\right) $. 
</s>