So, im trying to make use of this equation which im sure is well within my grasp but however I try to use it im doing it wrong, so was wondering if someone could offer a suggestion... Is just calculating the euclidian bit at the bottom I think im getting wrong (|| r - Ri ||). The summing over doesn't matter so much for the help I ask of you as the picture only describes a situation with 1 planet, but say I tried to apply that equation to that situation, trying to resolve for vx and vy - what would that look like ? For the record the best solution I have come up with is... v(x) = M (rx - R0) / (Sqrt(rx-R0))^3 I think, where x = y for the y component, but this creates a very erratic behaviour. As background this is being used in a eulers integral approximation method mixed with leapfrog, i.e. y(0) = 0 y(h) = y(0) + hF(0,v(x)) y(2h) = y(h) + hF(h, v(x)) thanks for any help