Hi everyone, I read an interesting tuto on mathematica. Pretty interesting for electrostatic: On the importance of being edgy in electrostatic I'm trying to find the trajectories of charged particles with random initial conditions in the vicinity of a uniformly charged disk with sharp edge. I defined the electrostatic field, then the gradient. After that I want to solve numerically the equations of the motion for several random conditions. I define 1st the potential as a sum of elliptic functions:

Q = 1; \[Epsilon] = 1; R = 1; r[x_, y_, z_] := Sqrt[(x)^2 + y^2]; 
p[x_, y_, z_] := Sqrt[(r[x, y, z] + R)^2 + (z)^2]; 
n[x_, y_, z_] := (4 r[x, y, z] R)/(r[x, y, z] + R)^2; 
m[x_, y_, z_] := (4 r[x, y, z] R)/p[x, y, z]^2; 
Ud[x_, y_, z_] := 
 2/p[x, y, 
   z] (p[x, y, z]^2*EllipticE[m[x, y, z]] - (r[x, y, z]^2 - R^2)*
     EllipticK[m[x, y, z]] - (z)^2*(r[x, y, z] - R)/(r[x, y, z] + R)
      EllipticPi[n[x, y, z], m[x, y, z]] ); 
Ui[x_, y_, z_] := 
 Piecewise[{{-2 \[Pi]*Sqrt[(z)*(z)], 
    r[x, y, z] \[LessSlantEqual] R}, {0, r[x, y, z] > R}}]; \[Phi][x_,
   y_, z_] := 
 1/(4 \[Pi] \[Epsilon])*Q/(\[Pi] R^2) (Ud[x, y, z] + Ui[x, y, z]);

Then I define the gradient:

Linedisk[{x_, y_, 
    z_}] := {-D[\[Phi]t[x, y, z], x], -D[\[Phi][x, y, z], 
     y], -D[\[Phi][x, y, z], z]};

After that I want to solve the equation of motion for several random conditions (here 5):

SeedRandom[111]; Show[{(*the charged disk *)
  Graphics3D[{Directive[GrayLevel[0.3], Specularity[Yellow, 25]], 
    Cylinder[{{0, 0, 0.01}, {0, 0, -0.01}}, 1]}],(*sample orbits *)
  Table[nds = 
    NDSolve[Join[
      Thread[{x''[t], y''[t], 
         z''[t]} == -Linedisk[{x[t], y[t], z[t]}]], 
      Thread[{x[0], y[0], z[0]} == 2 RandomReal[{-0.5, 0.5}/2, {3}]], 
      Thread[{x'[0], y'[0], z'[0]} == 
        RandomReal[{-0.5, 0.5}, {3}]]], {x, y, z}, {t, 0, 5}, 
     MaxSteps -> 50000]; 
   ParametricPlot3D[
    Evaluate[{x[t], y[t], z[t]} /. nds[[1]]], {t, 0.1, 
     nds[[1, 1, 2, 1, 1, 2]]}, 
    PlotStyle -> {ColorData["RedBlueTones"][RandomReal[]], 
      Tube[0.012]}, PlotRange -> All], {5}]}, PlotRange -> All]

After that I have many alarms like: Power::infy: "Infinite expression 1/0. encountered." or Infinity::indet: "Indeterminate expression ComplexInfinity+ComplexInfinity encountered." At the end, I have a 3D graph but I don't trust the results because of the alarms. There is a singularity at r=1 but normally it should be ok. I guess the problems are from the elliptic functions that are maybe too difficult and heavy for NDsolve. Any help will be appreciated to better understand the alarms from Mathematica. Thanks