How do I correctly introduce time into this equation? So, for the past few years it's been my goal to create an equation that would give me the position of an object in a gravitational field at time $t$, given it's initial position and velocity. At first the problem was that I didn't know enough to do the math. Now that I can do multivariable calculus I thought that problem would be solved, but I've just ended up running into a new problem. Please don't tell me how to solve it, but if you can give me a hint that would be great. Here's the set up for the problem:
A planet of mass $M$ (and radius = 0) is situated at the origin. I know that the magnitude of acceleration due to gravity is $$\frac{GM}{r^2}$$ so an object at $(x,y)$ will have acceleration $$a(x,y)= \frac{GM}{x^2+y^2},$$ or, as a vector, 
$$\overrightarrow{a}(x,y)= \left\langle \frac{GM}{x^2+y^2}\cos\theta, \frac{GM}{x^2+y^2}\sin\theta\right\rangle$$
$$= \left\langle \frac{GM}{x^2+y^2}\frac{x}{\sqrt{x^2+y^2}}, \frac{GM}{x^2+y^2}\frac{y}{\sqrt{x^2+y^2}}\right\rangle$$
$$= \left\langle \frac{GMx}{(x^2+y^2)^{3/2}}, \frac{GMy}{(x^2+y^2)^{3/2}}\right\rangle$$
So, here's where I'm stuck. I can integrate with respect to distance and get 
$$ W(x,y) = \left\langle -\frac{GM}{\sqrt{x^2+y^2}}, -\frac{GM}{\sqrt {x^2+y^2}}\right\rangle$$
which I think is a vector who's magnitude is the work done, but that doesn't tell me anything about time. I can integrate with respect to time, but that would give 
$$f(x,y)= \left\langle \frac{GMx}{(x^2+y^2)^{3/2}}t, \frac{GMy}{(x^2+y^2)^{3/2}}t\right\rangle$$
which... I mean is naïve at best. It doesn't take into account the change in position that happens over time. The only thing that I can think of to do is somehow find parametric equations where $x$ and $y$ are functions of $t$, but that's basically what I'm trying to do anyway.
Any ideas? I want to find an equation such that I can put in a location and velocity and the equation will tell me what path the object will take. Is that even possible?
 A: 
I can integrate with respect to time, but that would give
  $$f(x,y)= \left\langle \frac{GMx}{(x^2+y^2)^{3/2}}t, \frac{GMy}{(x^2+y^2)^{3/2}}t\right\rangle$$

No it doesn't because both $x$ and $y$ are functions of time, so you cannot use this "simple integration" method. What you need to do is go back to the differential equation,
$$
m\ddot{\mathbf r}=-\frac{GMm}{r^2}\hat{\mathbf r},\tag{1}
$$
and solve it for $r(t)$. Here we are using notation that means $\ddot{z}=d^2z/dt^2$ and $\dot z=dz/dt$ for some generic function $z$. When using vectors in non-Cartesian coordinates (e.g., cylindrical coordinates), it becomes a bit more complicated as the unit vectors themselves are time-varying. Solving (1) will require knowing some tricks (e.g., multiplying both sides by $\dot{r}$ and then integrating with respect to time), but it can be done.
I suggest solving this in one direction before trying it in two-dimensions first (e.g., in the $r$ direction before the $x,\,y$ directions). Solving it this way will show you the method (which is kinda complicated), doing it with the 2nd dimension will be a bit more straight-forward (since you have the experience)--but is also still pretty complicated.
Note also that right-hand side of (1) is the definition of $\mathbf g$, the gravitational acceleration:
$$
\mathbf g=-\frac{GM}{r^2}\hat{\mathbf r}.
$$
So (1) is really $m\mathbf a\equiv m\ddot{\mathbf r}\equiv\mathbf F=m\mathbf g$.
A: 
but if you can give me a hint that would be great.

Since the force depends on the radial distance only and points towards the origin, the angular momentum (assumed to point along the $z$ axis) is conserved.  This suggests that the appropriate coordinates are the spherical polar coordinates $(r, \phi)$ where
$$x = r \cos \phi$$
$$y = r \sin \phi$$
In these coordinates, the angular momentum is
$$l = mr^2\dot \phi = \mathrm{constant}$$
The radial equation of motion is
$$m\ddot r = mr\dot \phi^2 - \frac{GMm}{r^2}$$
Can you take it from here?
A: You have acceleration, which is a time derivative of velocity.  Integrating acceleration with respect to distance doesn't produce a work vector.  You can find the potential energy $U(x,y)=GMm/r$ which has the property that its gradient is the force field.  But the potential energy is a scalar, not a vector.
Your integral with respect to time is not correct, $x=x(t)$ and $y=y(t)$ are not constants in time, and you treated them like they were.
If you've studied multivariable calculus but haven't yet studied differential equations then study differential equations.
Also, start with motion in one dimension first (i.e. do a 1d problem), just because it is easier.  After you can do that easily, then go to higher dimensions.
A: Unfortunately, multivariable calculus isn't quite enough! Your goal is to solve a nonlinear system of ordinary differential equations, and that's a high goal.
Anyways, there are solutions to this. See an example of an implementation (disclaimer: It's mine) here: http://mathandcode.com/kepler/
The most useful references I found for it were the books Mathematical Methods of Classical Mechanics by Arnold, Mechanics by Landau and Lifshitz, and Classical Mechanics by Goldstein. Each of them have a solution of this problem - called "the Kepler problem". The solution outline goes something like this:


*

*Go from a differential equation in two variables $(x(t),y(t))$ to a differential equation in $r=\sqrt{x^2+y^2}$ only using the concept of an "Effective Potential".

*"Solve" this differential equation by abusing conservation of angular momentum (with which you know $\dot{\varphi}$), conservation of energy (with which you can find the quantity $\dot{r}^2+r^2 \dot{\varphi}^2$), and doing a clever change of variables so that you find $r$ as a function of $\varphi$ instead of as a function of time.

*Plug the equation for $r(\varphi)$ (which turns out to be the equation of an ellipse - of the form $\frac{p}{1+e \cos \varphi}$ for constants $p$ and $e$ called the parameter and eccentricity resp.) into the definition of angular momentum, which states $r^2 \dot{\varphi}$ is conserved and therefore gives you $\dot{\varphi}$

*Use the trick of separation of variables (which is the easiest ODE trick!) to find $\varphi(t)$ in terms of an integral. MAYBE then this integral can be evaluated in terms of special functions.


If this is all Greek to you, my advice is to pick up a book on mechanics (My favorites are Mathematical Methods of Classical Mechanics by Arnold, supplemented by The Theoretical Minimum [non-quantum] by Sussking and Hrabovsky).
