struct particle { vec3 p, v, f; }; std::vector<particle> boids; // Initialize N boids ... // ... // compute pairwise force for(int i=0; i<N; ++i) { for(int j=0; j<N; ++j) { if( i!=j ) { const vec3& pi = boids[i].p; const vec3& pj = boids[j].p; boids[i].f += force( norm(pi,pj)) / (pi-pj)/norm(pi-pj); } } } // integration for(int i=0; i<N; ++i) { boids[i].v = boids[i].v + dt * boids[i].f; boids[i].p = boids[i].p + dt * boids[i].v; }