Ways to decrease computational time?

Question

I have created a 2D axisymmetric simulation of the motion of charged particles in an electric field. The electric field is created by a needle attached to a top plate both of which are at -3800V with a grounded plate close below the needle. The needle is positioned above a hole in the plate which will allow particles to travel to the area below the grounded plate. The particles have an initial velocity, and in addition to the electric field, the forces acting on the particles include gravity, drag (approximated with Stokes Law), and the Coulomb force between particles.

It seems I cannot have more than two particles without it being too computationally intensive for my computer to solve (i.e. when there are 3 particles, it takes hours to solve for sol1 in the code below). Does anyone have any suggestions to improve the code and reduce computational time? Thanks!

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
q = -1.60217733*10^-19*10;(*particle charge*)
voltage = -3800;(*needle and top plate voltage*)
r1 = 0.0065; (*hole radius*)
r2 = 0.0365; (*domain radius*)
r3 = 0.00015; (*needle radius*)
z1 = 0.07; (*height of domain*)
z2 = 0.065; (*height of tip of needle*)
l = 0.005; (*distance between tip of needle and hole in middle plate*)
z3 = z2 - l; (*height of middle plate top surface*)
z4 = z3 - 0.0016; (*height of middle plate bottom surface*)
reg1 = ImplicitRegion[True, {{r, r3, r2}, {z, z2, z1}}]; (*region to the right of needle, above middle plate top surface*)
reg2 = ImplicitRegion[True, {{r, 0, r2}, {z, z3, z2}}]; (*region between tip of needle and middle plate top surface*)
reg3 = ImplicitRegion[True, {{r, 0, r1}, {z, z4, z3}}]; (*region to the left of middle plate*)
reg4 = ImplicitRegion[True, {{r, 0, r2}, {z, 0, z4}}]; (*region below middle plate bottom surface*)
region = RegionUnion[reg1, reg2, reg3, reg4]; (*merge all the previously created regions, space left out represents needle and middle plate*)
meshRefine[vertices_, area_] := area > 0.0000001;
mesh = ToElementMesh[DiscretizeRegion[region], MeshRefinementFunction -> meshRefine]; 
bc1 = {DirichletCondition[phi[r, z] == voltage, (z == z2 && 0 <= r <= r3)],
DirichletCondition[phi[r, z] == voltage, (z == z1 && r3 <= r <= r2)],
DirichletCondition[phi[r, z] == voltage, (r == r3 && z2 <= z <= z1)],
DirichletCondition[phi[r, z] == 0, (z == z3 && r1 <= r <= r2)],
DirichletCondition[phi[r, z] == 0, (z == z4 && r1 <= r <= r2)],
DirichletCondition[phi[r, z] == 0, (r == r1 && z4 <=z <=z3)]};  (*boundary conditions*)
sol = NDSolveValue[{1/r*D[r*D[phi[r, z], r], r] + 
      D[phi[r, z], z, z] == 0, bc1}, phi, {r, z} \[Element] mesh];
electricField = -Evaluate[Grad[sol[r, z], {r, z}]];
eforce = q*electricField;
dreal = 25*10^-9; (*real particle diameter in m*)
d = 7.5*10^-5; (*particle diameter in m (for simulation/visual purposes)*)
rho = 998; (*real particle density in kg/m^3*)
mu = 1.81*10^-5; (*air dynamic viscosity in Pa*s *)
volume = N[4/3*Pi*(dreal/2)^3]; (*real particle volume in m^3*)
mass = volume*rho;(*real particle mass in kg*)
VFR = 10/60/60*10^-6*0.001; (*volumetric flow rate through the needle in m^3/s*)
ivel = VFR/(Pi*r3^2); (*particle initial velocity out of the needle in m/s*)
numbodies = Round[r3/d]; (*number of particles that fit across the radius of the needle*)
gravity = -9.81; (*acceleration due to gravity in m/s^2*)
ke = 9*10^9; (*Coulomb constant*)
vel0 = Table[{0, -ivel}, numbodies];
pos0 = Riffle[Range[0, r3, r3/numbodies], z2]~Partition~2;
eqs = Table[{x[j]''[t] == 1/mass*Sum[((x[j][t] - x[i][t])/Sqrt[(x[j][t] - x[i][t])^2 
+ (y[j][t] - y[i][t])^2])*(ke*q^2/Sqrt[(x[j][t] - x[i][t])^2 + (y[j][t] - y[i][t])^2]),
 {i, Delete[Range[numbodies], j]}] + 1/mass*6.0*Pi*mu*dreal/2*- x[j]'[t]
 + 1/mass*eforce[[1]] /. {r -> x[j][t], z -> y[j][t]}, 
y[j]''[t] == gravity + 1/mass*Sum[((x[j][t] - x[i][t])/Sqrt[(x[j][t] - x[i][t])^2 
+ (y[j][t] - y[i][t])^2])*(ke*q^2/Sqrt[(x[j][t] - x[i][t])^2 + (y[j][t] - y[i][t])^2]), 
{i, Delete[Range[numbodies], j]}] + 1/mass*6.0*Pi*mu*dreal/2*-y[j]'[t] 
+ 1/mass*eforce[[2]] /. {r -> x[j][t], z -> y[j][t]}, x[j][0] == pos0[[j, 1]], 
y[j][0] == pos0[[j, 2]], x[j]'[0] == vel0[[j, 1]], y[j]'[0] == vel0[[j, 2]]}, {j, numbodies}];
vars = Flatten[Table[{x[j], y[j]}, {j, numbodies}]];
event = Table[{WhenEvent[x[j][t] == 0, {x[j]'[t] -> 0, y[j]'[t] -> 0}],
WhenEvent[x[j][t] == 0, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[x[j][t] == r2, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[y[j][t] == z1, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[y[j][t] == z3 && r1 <= x[j][t] <= r2, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[y[j][t] == z4 && r1 <= x[j][t] <= r2, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[x[j][t] == r1 && z4 <= y[j][t] <= z3, {x[j]'[t] -> 0, y[j]'[t] -> 0}]} /. j -> i, {i, numbodies}];
tfin = 25;
sol1 = NDSolve[{eqs, event}, vars, {t, 0, tfin}][[1]];
dp = DensityPlot[sol[r, z], {r, z} \[Element] mesh, ColorFunction -> "Rainbow", 
PlotLegends -> Automatic, PlotRange -> All, Frame -> False, AspectRatio -> Automatic];
frames = Table[Show[dp, ParametricPlot[Table[{x[j][t], y[j][t]} /. sol1, {j, numbodies}], {t, 0, tf}, 
PlotRange -> {{0, r2}, {0, z1}}, Axes -> False], Graphics[Table[{Hue[.35], 
Disk[{x[j][tf], y[j][tf]} /. sol1, d]}, {j, numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}];
ListAnimate[frames]