f = function(x,y) x-y oldpar = par(bg="black") plot(x,y,pch=".",col="white") for(i in 1:100) { #arrows(x,y,x+vx,y+vy,len=0.1) dx = outer(x,x,f) dy = outer(y,y,f) force = (dx*dx+dy*dy + 10)^(-1.5) diag(force) = 0; ax = apply(dx*force,2,sum) ay = apply(dy*force,2,sum) x = x+0.01*vx y = y+0.01*vy vx = vx + 0.01*ax vy = vy + 0.01*ay if(i%%10==0) plot(x,y,pch=".",col="white") } plot(x,y,pch=".",col="white")#,xlim=c(-10,10),ylim=c(-10,10),pch=".") par(oldpar)