NDSolve with arrays and Tables of Equations and WhenEvent

Question

I've been toying around with NDSolve for a while, and read through the website. By doing so I discovered that I could use it for vectors and arrays with much pleasure. So I wanted to write a simple simulator with multiple balls or particles interacting with Lennard Jones potentials. When I write this way I have no problem

partsim[p0_, v0_, mass_] := 
  Module[{rp},
  rp[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
  rp[t] /. NDSolve[
  Join[
   Table[
    {
     (*Initial COnditions*)
     rp[0][[i]] == p0[[i]],
     rp'[0][[i]]  == v0[[i]],

     (*Newton Equations*)  
     rp''[t][[i]] == (1/mass[[i]])*
       Sum[If[j != i, ljForce[i, j, rp[t]], 0], {j, 1, 
         nParticles}]
     (*Third newton law in the future*)

     }
    , {i, nParticles}](*End of Table*)  
   ](*End of Join*)
  , rp[t]
  , {t, 0, ftime}] // First
   ](*End of Function*);

The real problem arises when I'm trying to also add a WhenEvent for each of the particles present in the system. Since it's 3d I would like them to stay inside a bounding box and invert the sign of their velocity every time they touch the walls or at least come withing tolerance.

partsim[p0_, v0_, mass_] := 
  Module[{rp, vel},
  rp[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
  vel[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
  rp[t] /. NDSolve[
  Join[
   Table[
    {
     (*Initial COnditions*)
     rp[0][[i]] == p0[[i]],
     rp'[0][[i]] == vel[0][[i]] == v0[[i]],
     (*Newton Equations*)
     rp'[t][[i]] == vel[t][[i]],        
     rp''[t][[i]] == (1/mass[[i]])*
       Sum[If[j != i, pairForce[i, j, rp[t]], 0], {j, 1, 
         nParticles}]


     }
    , {i, nParticles}](*End of Table*)


   (*Discrete Events*)
   (*Bouncing on walls*)

   , Table[{
    WhenEvent[
     {Norm[rp[t][[u]] - wallcoord] < .01  , 
      Norm[rp[t][[u]]] <  .01}
     , vel[t][[u]] -> -vel[t][[u]]
     ](*End of WhenEvent*)
    }, {u, nParticles}]

   ](*End of Join*)

  , rp[t]
  , {t, 0, ftime}, DiscreteVariables -> {vel}] // First
   ](*End of Function*);

The error that appears seems to be related to the fact that u remains unevaluated, as attempting to Table it outside of the program will reveal. I have added DiscreteVariables -> {vel} because I thought that since we are subjecting velocity to a sudden change it would be necessary. Is there a way of adding a list of working WhenEvent, one for each particle, to the system of equations? Many thanks

Answer 1

In lieu of a working example in the current OP to work with, yes, this can be done. An example of two bouncing balls with differing velocity retention on bounce:

eventtab = 
 Table[With[{n = n}, 
   WhenEvent[y[n][t] == 0, y[n]'[t] -> (-0.95 + n/100) y[n]'[t]]], {n,
    2}]

NDSolve[{y[1]''[t] == -9.81, y[1][0] == 5, y[1]'[0] == 0, 
  y[2]''[t] == -4, y[2][0] == 3, y[2]'[0] == 0, eventtab}, {y[1], 
  y[2]}, {t, 0, 10}]

Plot[Evaluate[{y[1][t], y[2][t]} /. %], {t, 0, 10}]