Take a look at fvlib.  It's pretty darn fast.

Thanks Jeff. Have you already tried this library with 2d physics? There are only 2 examples with 3d models
and one with 3D point-spring..Can you share simple sketch with 2d physics?

My experience with it has only been with 3d, but here is what it might look like using 2d.

1. import volatileprototypes.fvlib.*;
2. Point[] psa;
3. Link[] lsa;
4. IntegratorVerlet vi;
5. BehaviorSpringRelaxation rs;
6. BehaviorConstantForce cf;
7. int steps=80;
8. int it=0;
9. float msec=0;
10. PFont font;
11. float strokew=3;
12. void setup() {
13.   ArrayList<Point> ps=new ArrayList();
14.   ArrayList<Link> ls=new ArrayList();
15.   size(640, 600);
16.   // Construct some random points.
17.   for (int i=0;i<2000;i++) {
18.     Point p=new Point(random(width)-width*.5, random(height)-height*.5, random(width)-width*.5, 1, .09, 200);
19.     ps.add(p);
20.   }
21.   Point p=new Point(random(width)-width*.5, random(height)-height*.5, random(width)-width*.5, 1, .09, 200);
22.   p.setU(true);
23.   ps.add(p);
24.   for (Point p1 : ps) {
25.     float nl=0;
26.     for (Point p2 : ps) {
27.       if (p1 != p2) {
28.         PVector pv = p2.get();
29.         pv.sub(p1);
30.         if (pv.mag() < 100 && nl < 6) {
31.           ls.add(new Link(p1, p2, .002));
32.           nl++;
33.         }
34.       }
35.     }
36.   }
37.   vi=new IntegratorVerlet(ps);
38.   rs=new BehaviorSpringRelaxation(ls).setFast(true);
39.   cf=new BehaviorConstantForce(ps, new PVector(0, .001, 0));
40.   // Convert ArrayLists to arrays for drawing performance.
41.   psa=new Point[ps.size()];
42.   lsa=new Link[ls.size()];
43.   ps.toArray(psa);
44.   ls.toArray(lsa);
45.   font = createFont("Arial", 12);
46.   for (int j=steps;j>=0;j--) {
47.     it++;
48.     msec*=.998;
49.     long start=System.currentTimeMillis(); //Benchmark
50.     rs.step();
51.     cf.step();
52.     vi.step();
53.     msec+=(float)(System.currentTimeMillis()-start)*.002; //Benchmark
54.   }
55. }
56. void draw() {
57.   float h=0;
58.   background(50, 100);
59.   pushMatrix();
60.   translate(width/2, height/2);
61.   Point p;
62.   fill(255);
63.   stroke(255, 30);
64.   strokeWeight(1);
65.   for (Link l : lsa) {
66.     line(l.p1.x, l.p1.y, l.p2.x, l.p2.y);
67.   }
68.   stroke(255, 70);
69.   strokeWeight(4);
70.   beginShape(POINTS);
71.   for (Point pt : psa) {
72.     vertex(pt.x, pt.y);
73.   }
74.   endShape();
75.   text("Points: "+ psa.length + " Springs: "+lsa.length, 10, height-49);
76.   text("Steps: "+ steps, 10, height-36);
77.   text("Total Loops: "+ steps*(psa.length+lsa.length), 10, height-23);
78.   text("Solver Time: "+ (int)(msec*steps) + "msec", 10, height-10);
79. }

Ok, but I can't figure out how to implement this on my falling boxes...

There was a huge bug with a fvlib module which has been fixed. So, if you download the latest version, the following code is possible:

Thanks Yiannis, now is more clearer..But it's possible with your library add initial velocity, parameters like density, friction and restitution (pbox2d), rotation, spin and boundaries??

To confirm, the elaborate collision geometry of Box2D is not geared to systems with such large numbers so going with a faster verlet physics engine (toxiclibs is another great choice) is a better option!

I fear that even with modern computers, you have to choose between the simplified and fast verlet integration engine (not much parameters) and the more sophisticated but much slower engine of JBox2D which can handle the parameters you wrote about...

I was thinking that it's possible to improve performance by putting boxes in sleeping status (no physics) when they reach the bottom..

Indeed as Daniel and Philho point out fvlib is focused almost exclusively in simulating particles and not solid bodies, at which Box2D excels, having much more features which fvlib does not of course aim to include. The roadmap of fvlib is right now all about improving performance and expanding particle-related behaviors further..

I can be wrong, but I believe that if you put boxes to sleep state, they become like the fixed part of the scene: the engine will no longer has to manage their movement, but still has to check for collisions with the still moving parts. Unless it has a super-smart algorithm eliminating checks against elements that are no longer reachable, or at least a good quadtree-like algorithm (checking only parts that are close together), it won't help much. Again, I can be wrong...
Now, perhaps one trick might consist in saying that if all the parts below a level are static, the program can take a shot of them, and replace all the parts with a big static block with an image of the eliminated parts...

In this sketch all boxes are initially in sleeping state. What I'm trying to do is add specific box to the box2d world when I click the mouse. All other boxes sleep till mouse is pressed. Please help...

download link

SOLVED !!!