Simulate "Newton's law of universal gravitation" using Box2D

It's pretty easy to implement:

for ( int i = 0; i < numBodies; i++ ) {

    b2Body* bi = bodies[i];
    b2Vec2 pi = bi->GetWorldCenter();
    float mi = bi->GetMass();

    for ( int k = i; k < numBodies; k++ ) {

        b2Body* bk = bodies[k];
        b2Vec2 pk = bk->GetWorldCenter();
        float mk = bk->GetMass();

        b2Vec2 delta = pk - pi;
        float r = delta.Length();
        float force = G * mi * mk / (r*r);

        delta.Normalize();
        bi->ApplyForce(  force * delta, pi );
        bk->ApplyForce( -force * delta, pk );
    }
}