Thank you myself!
I solved my problem by myself.
below is the thing I wanted to do. wierd motion of 2 charged particles!
Mathematica Code
r1 = {x1[t], y1[t], z1[t]};
r2 = {x2[t], y2[t], z2[t]};
v1 = Dt[r1, t];
v2 = Dt[r2, t];
deqns = Flatten[
{
Thread[
m1 Dt[v1, t] == (q1 q2)/(
4 Pi e0 ) ((r1 - r2)/(Norm[r1 - r2])^3 +
1/c^2 (v1\[Cross](v2\[Cross]((r1 - r2)/(Norm[r1 - r2])^3))))],
Thread[
m2 Dt[v2, t] == (q1 q2)/(
4 Pi e0 ) ((r2 - r1)/(Norm[r1 - r2])^3 +
1/c^2 (v2\[Cross](v1\[Cross]((
r2 - r1)/(Norm[r1 - r2])^3))))]
}
];
ics = {x1[0] == 0, y1[0] == 0, z1[0] == 0,
x1'[0] == 0, y1'[0] == 0, z1'[0] == 0,
x2[0] == 5, y2[0] == 0, z2[0] == 0,
x2'[0] == 0.1, y2'[0] == 0.1, z2'[0] == 0.1};
params = {q1 -> 10, q2 -> -10, m1 -> 1, m2 -> 1, e0 -> 1 , c -> 1};
soldp = First[NDSolve[{deqns, ics} /. params, {x1, y1, z1, x2, y2, z2}, {t, 0, 50}, MaxSteps -> Infinity]];
ParametricPlot3D[
Evaluate[{{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}} /. soldp],
{t, 0, 50}, PlotStyle -> {Red, Blue}, ImageSize -> Medium,
PlotLegends -> {"Trajectory1", "Trajectory2"}]
Animate[Graphics3D[{{PointSize[.025], {Red,
Point[{x1[t], y1[t], z1[t]}]}, {Blue,
Point[{x2[t], y2[t], z2[t]}]}} /. soldp,
{Red, Line[
Map[Function[Evaluate[{x1[#], y1[#], z1[#]} /. soldp]],
Range[0, t, .025]]]},
{Blue,
Line[Map[Function[Evaluate[{x2[#], y2[#], z2[#]} /. soldp]],
Range[0, t, .025]]]}
}, PlotRange -> {{0, 10}, {0, 10}, {0, 10}}, Axes -> True,
Ticks -> False, ImageSize -> 500], {t, 0, 50, .025},
SaveDefinitions -> True]