Considering a single particle moving more than one dimension (but still considering one particle), the equation of motions are:
F_i(\mathbf x) =m_i\,\ddot{x}_i
In general, the i^{th} component of the force \mathbf F depends from all the components, not only the i^{th} one, so there is not the direct relationship as in the one dimensional case (for which V would take the form \frac{\partial V_i}{\partial x_i}), but there is a potential V for which:
F_i(\mathbf x) = - \frac{\partial V(\mathbf x)}{\partial x_i}
This is the gradient which express the direction of the steepest descent. With this definition it is possible to prove conservation of energy for one particle, and it is done in the equivalent way:
E \equiv T + V = \sum \frac{1}{2}\,m\,\dot x ^2_i + V(\mathbf x)
Then taking the time derivative:
\begin{aligned} & \dot T = \sum_i \frac{\mathrm d}{\mathrm dt}\left(\frac{1}{2}\, m\,\dot x_i^2\right) = \sum_i \frac{m}{2}2\,\dot x_i\,\ddot x_i = \sum_i m\,\dot x_i\,\ddot x_i \\ & \dot V = \frac{\mathrm d}{\mathrm dt}V(x) = \sum_i \frac{\partial V}{\partial x_i} \dot x_i \\ & \dot E = \dot x_i \left( m\,\ddot x_i + \frac{\partial V}{\partial x_i} \right) = 0 \end{aligned}
In this case too using the definition of potential energy.
As example, it is possible to consider the following system which has the potential energy:
V = \frac{1}{2}k\left(x^2 + y^2 \right)
The first step is to change variables from x,y to r,\theta:
\begin{aligned} & x = r \cos(\theta) \\ & y = r \sin(\theta) \\ & x^2 + y^2 = r^2 \cos(\theta)^2 + r^2 \sin(\theta)^2 = r^2 \end{aligned}
Therefore:
V = \frac{1}{2}k\,r^2
It is possible to write the component of the force \mathbf F as \nabla \mathbf r in polar coordinates or in the two component F_x and F_y in Cartesian:
\begin{aligned} & \mathbf F = \nabla \mathbf r = \frac{\partial V}{\partial \mathbf r} = -k\,\mathbf r \\ & F_x = \frac{\partial V}{\partial x} = -k\,x \\ & F_y = \frac{\partial V}{\partial y} = -k\,y \end{aligned}
Then it is possible to write the equations of motion:
\begin{aligned} & F_x = m\, \ddot x = -k\,x \\ & \ddot x + \frac{k}{m}x = \ddot x + \omega^2 x = 0 \\ & F_y = m\, \ddot y = -k\,y \\ & \ddot y + \frac{k}{m}y = \ddot y + \omega^2 y = 0 \end{aligned}
These two second order equations can be written as a set of four first order equations:
\begin{aligned} & \dot x = u \\ & \dot y = v \\ & \dot u = -\omega^2 x \\ & \dot v = -\omega^2 y \end{aligned}
In matrix form [x_1,x_2,x_3,x_4]^T = [x,y,u,v]^T:
\begin{aligned} & \frac{\mathrm d}{\mathrm dt} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\omega^2 & 0 & 0 & 0 \\ 0 & -\omega^2 & 0 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} \\ & \dot{ \mathbf x} = \mathbf A \mathbf x \end{aligned}
To compute the solution it necessary to compute the characteristic polynomial as \det(\mathbf A -\lambda \, \mathbf I) =0:
\begin{aligned} & \begin{vmatrix} -\lambda & 0 & 1 & 0 \\ 0 & -\lambda & 0 & 1 \\ -\omega^2 & 0 & \lambda & 0 \\ 0 & -\omega^2 & 0 & \lambda \end{vmatrix} = 0 \\ & \lambda^4 + \omega^4 \end{aligned}
There are 2 repeated eigenvalue with multiplicity 2:
\begin{array}{ll} \lambda_1 = -i\,\omega & \mathbf v_1 = \begin{bmatrix} 0 & \dfrac{i}{\omega} & 0 & 1 \end{bmatrix}^T \\ \lambda_2 = -i\,\omega & \mathbf v_2 = \begin{bmatrix} \dfrac{i}{\omega} & 0 & 1 & 0 \end{bmatrix}^T \\ \lambda_4 = i\,\omega & \mathbf v_3 = \begin{bmatrix} 0 & - \dfrac{i}{\omega} & 0 & 1 \end{bmatrix}^T \\ \lambda_4 = i\,\omega & \mathbf v_4 = \begin{bmatrix} -\dfrac{i}{\omega} & 0 & 1 & 0 \end{bmatrix}^T \\ \end{array}
The generic solution can be written as:
\mathbf x(t) = \sum_{i=1}^4 c_i\mathbf v_ie^{\lambda_it}
Using the original vector:
\begin{bmatrix} x \\ y \\ u \\ v \end{bmatrix} = c_1e^{-i\omega t} \begin{bmatrix} 0 \\ \frac{i}{\omega} \\ 0 \\ 1 \end{bmatrix} + c_2e^{-i\omega t}\begin{bmatrix} \frac{i}{\omega} \\ 0 \\ 1 \\ 0 \end{bmatrix} + c_3 e^{i\omega t}\begin{bmatrix} 0 \\ - \frac{i}{\omega} \\ 0 \\ 1 \end{bmatrix} + c_4 e^{i\omega t} \begin{bmatrix} -\frac{i}{\omega} \\ 0 \\ 1 \\ 0 \end{bmatrix}
This gives:
\begin{aligned} & x(t) = \frac{i}{\omega} \left( c_2e^{-i\omega t} - c_4e^{i\omega t} \right) \\ & y(t) = \frac{i}{\omega} \left( c_1e^{-i\omega t} - c_3e^{i\omega t} \right) \\ & u(t) = c_2e^{-i\omega t} + c_4e^{i\omega t} \\ & v(t) = c_1e^{-i\omega t} + c_3e^{i\omega t} \end{aligned}
Using Euler formula:
\begin{aligned} & x(t) = \frac{i}{\omega} \left( c_2 \cos(-\omega\,t) + c_2 i\sin(-\omega\,t) - c_4 \cos(\omega\,t) - c_4 i\sin(\omega\,t) \right) \\ & y(t) = \frac{i}{\omega} \left( c_1 \cos(-\omega\,t) + c_1 i\sin(-\omega\,t) - c_3 \cos(\omega\,t) - c_3 i\sin(\omega\,t) \right)\\ & u(t) = c_2 \cos(-\omega\,t) + c_2 i\sin(-\omega\,t) + c_4 \cos(\omega\,t) + c_4 i\sin(\omega\,t) \\ & v(t) = c_1 \cos(-\omega\,t) + c_1 i\sin(-\omega\,t) + c_3 \cos(\omega\,t) + c_3 i\sin(\omega\,t) \end{aligned}
The constants can be determined with the initial conditions:
\begin{aligned} & x(0) = x_0 = \frac{i}{\omega} \left( c_2 \cos(0) + c_2 i\sin(0) - c_4 \cos(0) - c_4 i\sin(0) \right) = \frac{i}{\omega} (c_2 - c_4) \\ & y(0) = y_0 = \frac{i}{\omega} \left( c_1 \cos(0) + c_1 i\sin(0) - c_3 \cos(0) - c_3 i\sin(0) \right) = \frac{i}{\omega} (c_1 - c_3) \\ & u(0) = u_0 = c_2 \cos(0) + c_2 i\sin(0) + c_4 \cos(0) + c_4 i\sin(0) = c_2 + c_4\\ & v(0) = v_0 = c_1 \cos(0) + c_1 i\sin(0) + c_3 \cos(0) + c_3 i\sin(0) = c_1 + c_3 \end{aligned}
Solving:
\begin{aligned} & c_4 = u_0 -c_2 \\ & c_2 = \frac{1}{2} \left(u_0 + x_0 \frac{\omega}{i} \right) \\ & c_4 = \frac{1}{2} \left(u_0 - x_0 \frac{\omega}{i} \right) \\ & c_3 = v_0 -c_1 \\ & c_1 = \frac{1}{2} \left(v_0 + y_0 \frac{\omega}{i} \right) \\ & c_3 = \frac{1}{2} \left(v_0 - y_0 \frac{\omega}{i} \right) \end{aligned}
Simplifying as \cos(-x) = \cos(x) and \sin(-x) = -\sin(x):
\begin{aligned} & x(t) = \frac{i}{\omega} \left[ (c_2 - c_4) \cos(\omega\,t) - (c_2 + c_4) i\sin(\omega\,t) \right] \\ & y(t) = \frac{i}{\omega} \left[ (c_1 - c_3) \cos(\omega\,t) - (c_1 + c_3) i\sin(\omega\,t) \right]\\ & u(t) = (c_2 + c_4) \cos(\omega\,t) + (c_4 - c_2) i\sin(\omega\,t) \\ & v(t) = (c_1 + c_3) \cos(\omega\,t) + (c_3 - c_2) i\sin(\omega\,t) \end{aligned}
And inputting the constants:
\begin{aligned} & x(t) = \frac{i}{\omega} \left(x_0 \frac{\omega}{i} \cos(\omega\,t) - u_0 i\sin(-\omega\,t) \right) = x_0 \cos(\omega\,t) + \frac{u_0}{\omega} \sin(\omega\,t) \\ & y(t) = \frac{i}{\omega} \left[ (c_1 - c_3) \cos(\omega\,t) - (c_1 + c_3) i\sin(\omega\,t) \right] = y_0 \cos(\omega\,t) + \frac{v_0}{\omega} \sin(\omega\,t)\\ & u(t) = (c_2 + c_4) \cos(\omega\,t) + (c_4 - c_2) i\sin(\omega\,t) = u_0 \cos(\omega\,t) - x_0\omega \sin(\omega\,t) \\ & v(t) = (c_1 + c_3) \cos(\omega\,t) + (c_3 - c_2) i\sin(\omega\,t) = v_0 \cos(\omega\,t) - y_0 \omega \sin(\omega\,t) \end{aligned}
Therefore the equations of motions are:
\begin{array}{c} x(t) = x_0 \cos(\omega\,t) + \frac{u_0}{\omega} \sin(\omega\,t) \\ y(t) = y_0 \cos(\omega\,t) + \frac{v_0}{\omega} \sin(\omega\,t)\\ u(t) = u_0 \cos(\omega\,t) - x_0\omega \sin(\omega\,t) \\ v(t) = v_0 \cos(\omega\,t) - y_0 \omega \sin(\omega\,t) \end{array}
For this system is it possible to write the equations for the energy. The kinetic energy T is:
\begin{aligned} T & = \frac{1}{2}m \left( \dot x^2 + \dot y^2\right) = \frac{1}{2} m \left( \dot u^2 + \dot v^2\right) & = \frac{1}{2}m\left[\left( u_0 \cos(\omega\,t) - x_0\omega \sin(\omega\,t) \right)^2 + \left( v_0 \cos(\omega\,t) - y_0 \omega \sin(\omega\,t) \right)^2 \right] \end{aligned}
The potential energy V is:
\begin{aligned} V &= \frac{1}{2}k\left(x^2 + y^2 \right) \\ & = \frac{1}{2}k\left[\left( x_0 \cos(\omega\,t) + \frac{u_0}{\omega} \sin(\omega\,t) \right)^2 + \left( y_0 \cos(\omega\,t) + \frac{v_0}{\omega} \sin(\omega\,t) \right)^2 \right] \end{aligned}
The total energy is:
\begin{aligned} E = & T + V \\ = & \frac{1}{2}m \left[\left( u_0 \cos(\omega\,t) - x_0\omega \sin(\omega\,t) \right)^2 + \left( v_0 \cos(\omega\,t) - y_0 \omega \sin(\omega\,t) \right)^2 \right] \\ & + \frac{1}{2}k\left[\left( x_0 \cos(\omega\,t) + \frac{u_0}{\omega} \sin(\omega\,t) \right)^2 + \left( y_0 \cos(\omega\,t) + \frac{v_0}{\omega} \sin(\omega\,t) \right)^2 \right] \end{aligned}
It is possible to verify that the total energy is conserved:
\begin{aligned} \frac{\mathrm dE}{\mathrm dt} & = \frac{\mathrm dT}{\mathrm dt} + \frac{\mathrm dV}{\mathrm dt} = \frac{\mathrm d}{\mathrm dt} \left(\frac{1}{2}m\left( \dot u^2 + \dot v^2\right) + \frac{1}{2}k\left(x^2 + y^2 \right) \right) \\ & = \frac{1}{2}m\left( 2\dot x\ddot x + 2 \dot y \ddot y \right) + \frac{1}{2}k\left(2x \dot x + 2 y \dot y \right) = m\dot x\ddot x + m\dot y + \ddot y k x \dot x + ky \dot y \\ & = m\dot x\ddot x + m\dot y- m\dot x\ddot x - m\dot y = 0 \end{aligned}
as x = -\frac{m}{k}\ddot x and y = -\frac{m}{k}\ddot y.
The equation of motion can be rewritten as:
\begin{aligned} & x(t) = x_0 \cos(\omega\,t) + \frac{u_0}{\omega} \sin(\omega\,t) \\ & x(t) = A_x\cos(\omega\,t + \phi_x),\quad A_x = \sqrt{x_0^2 + \frac{u_0^2}{\omega^2}}, \quad \sin(\phi_x) = \frac{u_0}{A_x\omega} \\ & y(t) = y_0 \cos(\omega\,t) + \frac{v_0}{\omega} \sin(\omega\,t) \\ & y(t) = A_y \cos(\omega\,t + \phi_y), \quad A_y = \sqrt{y_0^2 + \frac{v_0^2}{\omega^2}}, \quad \sin(\phi_y) = \frac{v_0}{A_y\omega} \end{aligned}
Without loss of generality is always possible to translate such as \phi_x=0 and then \phi = \phi_y - \phi_x:
\begin{array}{c} x(t) = A_x\cos(\omega\,t) \\ y(t) = A_y\cos(\omega\,t + \phi) \end{array}
In order to obtain circlular orbits, there are conditions that need to be set on these quantities, as a circle requires the equation:
\begin{aligned} & x^2 + y^2 = r^2 \\ & A_x^2 \cos(\omega\,t)^2 + A_y^2 \cos(\omega\,t + \phi)^2 = A_x^2 \cos(\omega\,t)^2 + A_y^2 \sin(\omega\,t + \phi +\frac{\pi}{2})^2 = r^2 \\ & A_x^2 = A_y^2 = r^2 \\ & \phi = \frac{\pi}{2} \end{aligned}
Since these formula holds
\begin{aligned} & x_0 = A_x\,\cos(\phi_x) \quad \quad \frac{u_0}{\omega} = -A_x\,\sin(\phi_x) \\ & y_0 = A_y\,\cos(\phi_x) \quad \quad \frac{v_0}{\omega} = -A_x\,\sin(\phi_y) \end{aligned}
In the specific reference, with \phi_x = 0 and \phi_y = \frac{\pi}{2}:
\begin{array}{cc} x_0 = A_x\,\cos(0) & \dfrac{u_0}{\omega} = -A_x\,\sin(0) \\ x_0 = A_x & u_0 =0\\ y_0 = A_y\,\cos\left(\dfrac{\pi}{2}\right) & \dfrac{v_0}{\omega} = -A_y\,\left(\dfrac{\pi}{2}\right) \\ y_0 = 0 & v_0 = A_y \end{array}
In case that A_x \neq A_y the orbits are ellipsis; if \phi \neq \frac{\pi}{2} the axis of the ellipsis is not aligned with the major axes with the x and y axes. This is the isotropic case in which the potential is proportional to \mathbf r (or equivalently that k_x = k_y or \omega_x = \omega_y).
There is also another case, the anisotropic one in which the potential is not radially constant:
V = \frac{1}{2}k_x x^2 + k_y y^2
The equations of motion are similar to the previous case:
\begin{aligned} & F_x = m\, \ddot x = -k_x\,x \\ & \ddot x + \frac{k}{m}x = \ddot x + \omega_x^2 x = 0 \\ & F_y = m\, \ddot y = -k_y\,y \\ & \ddot y + \frac{k}{m}y = \ddot y + \omega_y^2 y = 0 \end{aligned}
To arrive to the solution there is a similar system:
\begin{array}{c} x(t) = x_0 \cos(\omega_x\,t) + \frac{u_0}{\omega_x} \sin(\omega_x\,t) \\ y(t) = y_0 \cos(\omega_y\,t) + \frac{v_0}{\omega_y} \sin(\omega_y\,t)\\ u(t) = u_0 \cos(\omega_x\,t) - x_0\omega_x \sin(\omega_x\,t) \\ v(t) = v_0 \cos(\omega_y\,t) - y_0 \omega_y \sin(\omega_y\,t) \end{array}
or equivalently:
\begin{array}{c} x(t) = A_x\cos(\omega_x\,t) \\ y(t) = A_y\cos(\omega_y\,t + \phi) \end{array}
There is the convention that the time start at the turning point of x and \phi = \phi_y - \phi_x. In the anisotropic case the shape of the orbits cannot easily identified and depends from the frequencies and the initial conditions.