public double getF (int id,double x[],double y[]) {
  for(int i=0;i<n;i++){
    dx=x[id]-x[i]; dy=y[i]-y[id];
//    if(r>R3)f2=0.;
//    else if(r>R) f2=F0*(1.5-0.5*r/R)/r;
    else f2=F0+(R-r)*(Fmax-F0)/R;
//     f2=F0*R/r;
//    if(r     f+=f2*dx/r;
  return f/n;

I have tested the above code with different models. (some lines has been removed, some are just comment out!)
As far as I can remembered: I did tried with Lennard-Jones potential , but I did not get a good result. (did not find the good parameter space values for the simulation).
You can tried with different model by yourself. This is the beauty of using EJS to build and test you own model.

The above code assume F(r)=-Fo/r; in the vector [b]r[/b] direction.
Then, the x component is calculated as Fx(r)=F(r)*x/r=-Fo*x/r[sup]2[/sup]; (similar for Fy(r)=F(r)*y/r;)
Actually, this is  Columb's Law in 2D.
For Columb's Law in 3D: The formula is the famous F(r)=-F0/r[sup]2[/sup];

You can understand the above relations with [u]Gauss Law[/u]. The model assume space is uniform.
In 3D: The surface area at r is 4*?*r[sup]2[/sup] so  force is proportional to  r[sup]-2[/sup] (or 1/r[sup]2[/sup]).
In 2D: The surface area at r is 2*?*r so  force is proportional to  r[sup]-1[/sup] (or 1/r).
Just like the B field at distance r from a straight current line is proportional to 1/r;