There's no way I can reproduce the run, since you haven't posted full code; I'm doing this from visual inspection.
Among other things, you've left variables k, T, m undefined, so I have no idea of the magnitude of the velocities. Is the upper limit high enough for a particle to skip totally past a collision? If v*dt can be > tol*2, you can totally miss collisions. When I code things like this, I make sure that tol is bounded by the maximum change.
Other than that, your position update is correct.
What happens when you run this with a single particle? Can you get it to bounce off all three walls? What happened when you hand-coded two particles headed straight at one another along a single dimension? I would like to think that you tested these basic functions before you loaded in 10000 particles.
I'm not sure what L is in your code; the only rational idea I have is that it's the upper limit of the box dimensions, which I'm reading as 1.2^10-6 on your plot.
The biggest problem that I see is that your only collision check is against the x+ side of the box. You haven't check the other dimensions or the lower limit (0?) for anything, and there's no comparison of one particle to another. As a result, the only change you'll get is to have roughly half your particles hit the right wall (x+) and reverse direction. Other than that, everything will drift forever, with ever-decreasing x values.
HANDLE A BOUNCE
First of all, get your velocity and position expressions in synch with one another: each of these should be a triple (tuple or list), so you can index them together. You currently have them as the second index of your particle list (good idea), but three separate variables elsewhere (poor idea). Instead, make newVel a list, so you can simply loop through the particles and update:
for dim in range(3):
if r[i, dim] < tol or # hit the lower wall
r[i, dim] > b - tol: # hit the upper wall
v[i, dim] = -v[i, dim]
Also, please just update the positions appropriately; no need to involve the local, temporary variables:
for dim in range(3):
r[i, dim] += v[i, dim] * dt
WRITE SOME SERVICE FUNCTIONS
It often helps to write some general-use functions, just to get them out of the way when you're writing your main code. At the moment, you're handling a lot of details, when you should be worrying about only one technique at a time. For instance, since you're going to be handling particle collisions, you'll want a distance computation. Don't keep this in the middle of your active code; just write a function to do it. For instance:
def dist(a, b): # a, b are positions: triples of coordinates
return math.sqrt(sum([(a[dim] - b[dim])**2 for dim in range(3)]))
I won't comment on how to implement the particle collisions, because you haven't shown any effort to solve the problem yet, and you're not yet at a point to attack that extension. First, get this part working: one particle bouncing around, then two particles bouncing around, then maybe 4 particles to show that you can do it for an arbitrary quantity of particles.
Once you have that, you can worry about particle collisions. When you get to that point, make a good attempt; if you get stuck, post a new question.
