3D Parametric Plotting for N Charged Particles in a Nonuniform Electric Field

I have a code (very bottom of post) which plots a static electric field as it passed through some metal sheet with an aperture in it.

I seek to observe the motion of some N charged particles (point charges) randomly dispersed near the aperture-field (seen above). I have created a set of equations for the N particles in x[t], y[t] and z[t] as they sit in the electric field, and can solve the set of differential equations which exist for each particle.

When it comes time to plot the motion of the particles and the electric field/aperture (seen above), I am running into some problems in the last few lines of the code

frames = Table[ Show[v, ParametricPlot3D[ Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j, numbodies}], {t, 0, tf}, PlotRange -> Automatic, Axes -> False], Graphics3D[ Table[{Hue[.35], Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.0005]}, {j, numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}]; ListAnimate[frames] 

which produces these errors:

Here is the full code:

Clear["Global`*"]; Needs["NDSolve`FEM`"] q = 1.60217733*10^-19*10;(*Net ion Charge*) R = Import["https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=1"]; RegionDimension[R]; M = BoundaryMeshRegion[MeshCoordinates[R], MeshCells[R, 2]]; RegionDimension[M]; Volume[M]; r = RegionDifference[ RegionDifference[ RegionDifference[Cuboid[{0, 0, -0.5}, {2, 2, 0.5}], M], Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]], Cuboid[{0, 0, -0.5}, {2, 2, -0.4}]]; ToElementMesh[r]["Wireframe"]; pol = -1; V0 = 4000; sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0, DirichletCondition[ V[x, y, z] == -pol* V0/2, (0.4 <= z <= 0.5) && (0 <= y <= 2) && (0 <= x <= 2)], DirichletCondition[ V[x, y, z] == pol*V0/2, (0.0071 <= z <= 0.0072) && (0 <= y <= 2) && (0 <= x <= 2)], DirichletCondition[ V[x, y, z] == 0, (0 <= z <= 0.0070) && (0 <= y <= 2) && (0 <= x <= 2)], DirichletCondition[ V[x, y, z] == 0, (-0.5 <= z <= -0.4) && (0 <= y <= 2) && (0 <= x <= 2)]}, V, {x, y, z} \[Element] r] electricField[x_, y_, z_] = -Grad[sol[x, y, z], {x, y, z}]; v = Show[VectorPlot3D[ electricField[x, y, z], {x, 0.5, 1}, {y, 0.5, 1}, {z, -0.25, -0.00001}, PlotTheme -> "Detailed", ColorFunction -> "Rainbow", PerformanceGoal -> "Quality", VectorScale -> 0.05, VectorPoints -> 7], M]; eforce1 = q*electricField1[x, y, z]; vecForce = Show[VectorPlot3D[ q*electricField[x, y, z], {x, 0.5, 1}, {y, 0.5, 1}, {z, -0.25, -0.00001}, PlotTheme -> "Detailed", ColorFunction -> "Rainbow", PerformanceGoal -> "Quality", VectorScale -> 0.05, VectorPoints -> 7], M] mass = 6.52*10^-8;(*particle mass in kg/m^3*) numbodies = 3; vel0 = Table[ Partition[{RandomReal[{-0.0001, 0.0001}], RandomReal[{-0.0001, 0.0001}], RandomReal[{-0.0001, 0.0001}]}, 1], numbodies] pos0 = Table[ Partition[{RandomReal[{0, 2}], RandomReal[{0, 2}], RandomReal[{0.00001, 1}]}, 1], numbodies] eqs = Table[{x1[j]''[t] == 1/mass*eforce1[[1]] /. {x -> x1[j][t]}, y1[j]''[t] == 1/mass*eforce1[[1]] /. {y -> y1[j][t]}, z1[j]''[t] == 1/mass*eforce1[[1]] /. {z -> z1[j][t]}, x1[j][0] == pos0[[j, 1]], y1[j][0] == pos0[[j, 2]], z1[j][0] == pos0[[j, 3]], x1[j]'[0] == vel0[[j, 1]], y1[j]'[0] == vel0[[j, 2]], z1[j]'[0] == vel0[[j, 3]]}, {j, numbodies}]; vars = Flatten[Table[{x1[j], y1[j], z1[j]}, {j, numbodies}]] event = Table[{WhenEvent[ x1[j][t] == 0, {x1[j]'[t] -> 0, y1[j]'[t] -> 0, z1[j]'[t] -> 0}], WhenEvent[ y1[j][t] == 0, {x1[j]'[t] -> 0, y1[j]'[t] -> 0, z1[j]'[t] -> 0}], WhenEvent[ z1[j][t] == 0, {x1[j]'[t] -> 0, y1[j]'[t] -> 0, z1[j]'[t] -> 0}]} /. j -> i, {i, numbodies}]; tfin = 15 sol1 = NDSolve[{eqs, event}, vars, {t, 0, tfin}][[1]] plotXZ = ContourPlot[sol[x, 0.75, z], {x, 0, 2}, {z, -0.4, 0.1}, ContourShading -> Automatic, ColorFunction -> "Rainbow", Contours -> 100]; frames = Table[ Show[v, ParametricPlot3D[ Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j, numbodies}], {t, 0, tf}, PlotRange -> Automatic, Axes -> False], Graphics3D[ Table[{Hue[.35], Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.0005]}, {j, numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}]; ListAnimate[frames] ```