Oribit integrator for a logarithmic potential

Hello! Right know I'm trying to make an orbit integrator for solving a logarithmic potential with the form:
\begin{equation}
\Phi= \frac{v_0^2}{2} ln(x^2+ \frac{y^2}{u^2} + r_0^2)
\end{equation}
where v0, u, and r0 are constants
My approach is to use,
\begin{equation}
\ddot{q} = -\bigtriangledown \Phi
\end{equation}
Then the system equations,
\begin{equation}
\ddot{x} = -v_o^2 \frac{x}{x^2+ \frac{y^2}{u^2} + r_0^2}
\end{equation}
\begin{equation}
\ddot{y} = -\frac{v_o^2}{u^2} \frac{y}{x^2+ \frac{y^2}{u^2} + r_0^2}
\end{equation}
My guess is that in order to solve for x and y using Runge Kutta or leapfrog, I need to decouple the system, but I don't know how to do so.

If by decouple the system, you mean separate the variables, it probably is not possible. But Runge-Kutta can be used to integrate the equations of motion as is. The Hamiltonian formulation (four equations with first derivatives) is usually easier.

By doing the Hamiltonian approach I still get equations (3) and (4) above, and the other two are are apparently of no use.
The problem is that I don't know how (if possible) to adapt the Runge-Kutta using two dependent variables (x,y) and the independent one (t).