For two bodies this is relatively easy as the equations of motion describe a conic (an ellipse for a closed orbit, a hyperbola for an "open" orbit). You can use the vis viva equation to get the parameters of the orbit (semi major axis etc) from the given initial conditions, and the rest follows.
For an ellipse, you can also express the position as a function of time as
$$v_x = a \sin \omega t\\
v_y = b \cos \omega t$$
where your challenge is to find $a$, $b$ and $\omega$. For this you would like to know the velocity and distance at the furthest point, but you can in principle solve this for any place along the orbit.
For example, for a satellite in elliptical orbit above a massive planet (mass $M$), for given initial velocity $v$ in the radial direction, at a distance $r$ from the center of the planet, we can know the following (see http://en.wikipedia.org/wiki/Vis-viva_equation):
$$v^2 = GM\left(\frac{2}{r}-\frac{1}{a}\right)$$
In this expression, $a$ is the semimajor axis. We can calculate it from the given initial conditions,
$$\frac{1}{a} = \frac{2}{r}-\frac{v^2}{2GM}$$
Now we need $b$; again, following the wikipedia page quoted above, and for the conditions given (v perpendicular to radial vector at distance r, so we know angular momentum of the orbit $L=mvr$), we can write
$$b = vr\sqrt{\frac{a}{GM}}$$
With $a$ and $b$ calculated, we just need to find $\omega$ which follows again from the initial conditions, since we know
$$\omega = \frac{v}{r}$$
at the initial position. And with $a$, $b$ and $\omega$ known you can substitute into the above equations.
Incidentally it is more practical, when you have more than one gravitational object, to do the integration "properly". I recommend using a fourth-order Runge-Kutta integration, which is quite stable for this kind of problem. You can see http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf for some ideas on how to do this - it's a bit hard to follow. I actually wrote a piece of code that did this in Visual Basic a little while ago - it was actually taking into account both gravity and atmospheric drag. I reproduce it here without warranty...
Option Explicit
Const gm = 9.8 * 6300000# * 6300000# ' approximate value of GM on earth surface
' the # signifies a double precision "integer"
Const C_d = 0.5 * 0.03 * 0.03 * 0.2 ' drag coefficient for particular object
Const R_e = 6300000# ' approximate radius of earth
Const mass = 25 ' mass of object
Const k1 = 12000# ' constants used to fit density of atmosphere
Const k2 = 22000# ' again, approximate!
Function rho(h)
' return density as a function of height
If h < 0 Then rho = 1.2 Else rho = 1.225 * Exp(-(h / k1 + (h / k2) ^ (3 / 2)))
End Function
Sub accel(x, y, vx, vy, ByRef ax, ByRef ay)
' compute acceleration of object at location x,y with velocity vx, vy
' return values in ax, ay
Dim r As Double, v2 As Double, v As Double, r3 As Double, h As Double
Dim density
r = Sqr(x * x + y * y) ' sqrt() in most languages... this is VBA
v2 = vx * vx + vy * vy
v = Sqr(v2)
r3 = 1# / r ^ 3
density = rho(r - R_e)
' don't need the second term if you don't care about drag!
ax = -gm * x * r3 - vx ^ 3 * C_d * density / (v * mass)
ay = -gm * y * r3 - vy ^ 3 * C_d * density / (v * mass)
End Sub
Function rk(ByRef x, ByRef y, ByRef vx, ByRef vy, dt)
' implement one Runge-Kutta fourth order stop
' this requires calculating the acceleration at verious locations
' for a single "step" in the algorithm
Dim ax As Double, ay As Double
Dim kx1, kx2, kx3, kx4, ky1, ky2, ky3, ky4, lx1, lx2, lx3, lx4, ly1, ly2, ly3, ly4, dt2
' get acceleration at initial point
accel x, y, vx, vy, ax, ay
' half time step:
dt2 = dt / 2
kx1 = dt2 * ax
ky1 = dt2 * ay
lx1 = dt2 * vx
ly1 = dt2 * vy
' get acceleration at new location
accel x + lx1, y + ly1, vx + kx1, vy + ky1, ax, ay
kx2 = dt2 * ax
ky2 = dt2 * ay
lx2 = dt2 * (vx + kx1)
ly2 = dt2 * (vy + ky1)
' get acceleration at half way point:
accel x + lx2, y + ly2, vx + kx2, vy + ky2, ax, ay
kx3 = dt * ax
ky3 = dt * ay
lx3 = dt * (vx + kx2)
ly3 = dt * (vy + ky2)
' compute acceleration for third combination of position and velocity:
accel x + lx3, y + ly3, vx + kx3, vy + ky3, ax, ay
kx4 = dt2 * ax
ky4 = dt2 * ay
lx4 = dt2 * (vx + kx3)
ly4 = dt2 * (vy + ky3)
' compute best approximation to new x, y, vx, vy
x = x + (lx1 + 2# * lx2 + lx3 + lx4) / 3#
y = y + (ly1 + 2# * ly2 + ly3 + ly4) / 3#
vx = vx + (kx1 + 2# * kx2 + kx3 + kx4) / 3#
vy = vy + (ky1 + 2# * ky2 + ky3 + ky4) / 3#
End Function
Calling the RK function multiple times you can now track the orbit, and you will find it is quite stable. You can easily add multiple massive objects (just update the acceleration formula) and it will still work (although orbits will become chaotic).
