Expressing position due to gravitational acceleration as a 3-Dimensional differential equation I know that the force of gravity is $F = \dfrac{G  m_1   m_2}{r^2}$. Now assume in a one-dimensional system there are two masses in the universe, a planet and an object. The object would have an acceleration of $A = \dfrac{G  M}{r^2}$, assuming that m is the mass of the planet. The position of the object can then be graphed with the differential equation $\frac{d^2x}{dt^2}=-\frac{GM}{x^2}$, with x being position and t being time, and negative as the object falls towards the planet. This differential equation can be solved with initial position and velocity values.



My question is now how to transform this to 2 (and 3) dimensions. Assume the same situation but now the object has an x and y coordinate. I can say that acceleration in the x direction equals $A_x = \dfrac{G  M}{r^2} \dfrac{x}{r}$, with r being the radius ($\sqrt{y^2+x^2}$. I can also say that acceleration in the y direction equals $A_y = \dfrac{G  M}{r^2} \dfrac{y}{r}$.



These equations work when I graph them, but because they rely on each other (the x position is needed for y-acceleration and vice versa), I'm unable to find a way to solve the differential equations and find x and y as functions.
Does anyone have any ideas? If I need to clarify myself better, please ask. Thank you for your time!
 A: 
My question is now how to transform this to 2 (and 3) dimensions

I am not sure where your problem is in doing this. In 3 dimensions using
$$\mathbf{x}_1 = \begin{pmatrix}x_1\\y_1\\z_1\end{pmatrix}, \mathbf{x}_2 = \begin{pmatrix}x_2\\y_2\\z_2\end{pmatrix}$$
to denote the position of $m_1$ and $m_2$ respectively, the equations of motion are:
$$m_1\ddot{\mathbf{x}}_1=\mathbf{F}_{21}\\
m_2\ddot{\mathbf{x}}_2=\mathbf{F}_{12}$$
where $\mathbf{F}_{21}$ is the gravitational force of $m_2$ on $m_1$ and $\mathbf{F}_{12}$ is the gravitational force of $m_1$ on $m_2$. Because of Newton's third law, $\mathbf{F}_{21}=-\mathbf{F}_{12}$.
This is a somewhat complicated system of non-linear differential equations. However if you think about the symmetries of the system it should be clear that the motion should not depend on the absolute position of the masses in space but only on their relative position. Mathematically this simplification can be achieved by adding/subtracting the equations from one another. Adding the equations of motions gives:
$$0=m_1\ddot{\mathbf{x}}_1+m_2\ddot{\mathbf{x}}_2= (m_1+m_2)\ddot{\mathbf{R}}$$, 
where $\mathbf{R}\equiv (m_1\mathbf{x}_1+m_2\mathbf{x}_2)/(m_1+m_2)$ is the center of mass coordinate. This equation expresses the conservation of momentum. 
In the same way for the difference of the coordinates you get the equation of motion:
$$\ddot{\mathbf{r}}\equiv \ddot{\mathbf{x}}_2 - \ddot{\mathbf{x}}_1 = \left(\frac{1}{m_1}+\frac{1}{m_2}\right)\mathbf{F}_{12}$$ which can be rewritten with the reduced mass, $\mu=m_1m_2/(m_1+m_2)$ as
$$\mu\ddot{\mathbf{r}}=\mathbf{F}_{12}=-\frac{Gm_1m_2}{r^2}\hat{\mathbf{r}}$$
where $r=|\mathbf{r}|$ is the distance between the masses and $\hat{\mathbf{r}}=\mathbf{r}/r$ is the unit vector in the direction of $\mathbf{r}$. This equation is the starting point for the discussion/calculation of the motion.
Outline of solution
Finding a general analytical solution for $\mathbf{r}(t)$ is difficult (impossibe or at least ugly requiring special functions?). Still you can get many important results without this solution. 
Using symmetry or angular momentum conservation you can argue that the motion will be in a plane (set by the initial velocity vectors), so that the problem simplifies to a two dimensional problem.
Since the force is radial it makes sense to solve the problem in polar coordinate:
$$\mathbf{r}=r\begin{pmatrix}\cos\theta\\ \sin\theta\end{pmatrix}$$
Even then solving the problem is hard work. The classical method is to substitute $u=1/r$ and to look for solutions for $u(\theta)$. You might want to look up Binet equation.
A: The easiest thing to do is to start from the gravitational potential
$$
\phi=-\frac{Gm_1m_2}{\sqrt{x^2+y^2+z^2}}
$$
and use $\vec F=-\vec\nabla \phi$, or in components:
$$
F_x=-\frac{\partial \phi}{\partial x}= 
-\frac{Gm_1m_2\,x}{(x^2+y^2+z^2)^{3/2}} = m_2 A_x = m_2 \frac{d^2 x}{dt^2}
$$
with a similar approach for your other components.
