(* Relativistic Gravity & Charge Simulator. Yukterez, Version 1.8.2014 *)

(* Sonne *)
M1=mSol;
Q1=0; 
P1x=+4.767411813724087*^+5 km;
P1y=-4.881392129931528*^+4 km;
P1z=-2.184128405553382*^+4 km;
v1x=+5.632851183659290*^-3 km/sek;
v1y=+1.020293715934504*^-2 km/sek;
v1z=-1.481762989445804*^-4 km/sek;

(* Merkur *)
M2= 3.302*^23 kg;
Q2=0;
P2x=+5.357892741028007*^+7 km;
P2y=-1.147221840127503*^+6 km;
P2z=-4.983469907397524*^+6 km;
v2x=-8.396224055586826*^+0 km/sek;
v2y=+5.089273470880620*^+1 km/sek;
v2z=+4.928244620166743*^+0 km/sek;

(* Venus *)
M3=4.8685*^24 kg;
Q3=0;
P3x=-3.828375486384919*^+7 km;
P3y=+1.001776253505679*^+8 km;
P3z=+3.588675833902176*^+6 km;
v3x=-3.277640635284602*^+1 km/sek;
v3y=-1.281675204600856*^+1 km/sek;
v3z=+1.715864558132450*^+0 km/sek;

(* Erde *)
M4= 5.97219*^24kg; 
Q4=0;
P4x=-1.436305243395423*^+8 km;
P4y=-4.046527983363845*^+7 km;
P4z=-2.022577828499544*^+4 km;
v4x=+7.558464309586846*^+0 km/sek;
v4y=-2.877403274780793*^+1 km/sek;
v4z=-2.203233704649544*^-4 km/sek;

(* Schweif *)
Tmin=yr/12;

(* Simulationszeit in sek *)
τ=3yr;

(* Eckdaten *)
G=667384/10^16(m^3)(kg^(-1))(sek^(-2)); (*Gravitationskonstante*) 
ε0=8.854187817*^-12;  (* Feldkonstante *) 
Au=149597870690m; (*Astronomische Einheit*) 
mSol=1.988544*^30 kg; (*Sonnenmasse*) 
c=299792458m/sek; (* Licht *)
kg=1; m=1; sek=1;km=1000m;  yr=31536000sek;

(* Relativistische Masse, longitudinal und transversal *)

γ1x:=((((Sqrt[1-(((vx1[t])^2)/(c^2))])))^3)(Sqrt[1-((((vy1[t])^2)+((vz1[t])^2))/(c^2))]); γ2x:=((((Sqrt[1-(((vx2[t])^2)/(c^2))])))^3)(Sqrt[1-((((vy2[t])^2)+((vz2[t])^2))/(c^2))]);γ3x:=((((Sqrt[1-(((vx3[t])^2)/(c^2))])))^3)(Sqrt[1-((((vy3[t])^2)+((vz3[t])^2))/(c^2))]); γ4x:=((((Sqrt[1-(((vx4[t])^2)/(c^2))])))^3)(Sqrt[1-((((vy4[t])^2)+((vz4[t])^2))/(c^2))]);

γ1y:=(Sqrt[1-((((vx1[t])^2)+((vz1[t])^2))/(c^2))])((((Sqrt[1-(((vy1[t])^2)/(c^2))])))^3); γ2y:=(Sqrt[1-((((vx2[t])^2)+((vz2[t])^2))/(c^2))])((((Sqrt[1-(((vy2[t])^2)/(c^2))])))^3);γ3y:=(Sqrt[1-((((vx3[t])^2)+((vz3[t])^2))/(c^2))])((((Sqrt[1-(((vy3[t])^2)/(c^2))])))^3); γ4y:=(Sqrt[1-((((vx4[t])^2)+((vz4[t])^2))/(c^2))])((((Sqrt[1-(((vy4[t])^2)/(c^2))])))^3);

γ1z:=(Sqrt[1-((((vx1[t])^2)+((vy1[t])^2))/(c^2))])((((Sqrt[1-(((vz1[t])^2)/(c^2))])))^3); γ2z:=(Sqrt[1-((((vx2[t])^2)+((vy2[t])^2))/(c^2))])((((Sqrt[1-(((vz2[t])^2)/(c^2))])))^3);γ3z:=(Sqrt[1-((((vx3[t])^2)+((vy3[t])^2))/(c^2))])((((Sqrt[1-(((vz3[t])^2)/(c^2))])))^3); γ4z:=(Sqrt[1-((((vx4[t])^2)+((vy4[t])^2))/(c^2))])((((Sqrt[1-(((vz4[t])^2)/(c^2))])))^3);


(* Kraft *)

Fx1:=(((G  (m1/(γ1x)) (m2/(γ2x)) (x2[t]-x1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((G  (m1/(γ1x)) (m3/(γ3x)) (x3[t]-x1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((G  (m1/(γ1x)) (m4/(γ4x)) (x4[t]-x1[t]))/(Sqrt[((((((x4[t]-x1[t]))^2)+(((y4[t]-y1[t]))^2)+(((z4[t]-z1[t]))^2)))^3)]))+

((-q1*q2/(4π ε0 )(x2[t]-x1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((-q1*q3/(4π ε0 )(x3[t]-x1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((-q1*q4/(4π ε0 ) (x4[t]-x1[t]))/(Sqrt[((((((x4[t]-x1[t]))^2)+(((y4[t]-y1[t]))^2)+(((z4[t]-z1[t]))^2)))^3)])));


Fy1:=(((G  (m1/(γ1y)) (m2/(γ2y)) (y2[t]-y1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((G (m1/(γ1y)) (m3/(γ3y)) (y3[t]-y1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((G  (m1/(γ1y)) (m4/(γ4y)) (y4[t]-y1[t]))/(Sqrt[((((((x4[t]-x1[t]))^2)+(((y4[t]-y1[t]))^2)+(((z4[t]-z1[t]))^2)))^3)]))+

((-q1*q2/(4π ε0 )(y2[t]-y1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((-q1*q3/(4π ε0 )(y3[t]-y1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((-q1*q4/(4π ε0 )(y4[t]-y1[t]))/(Sqrt[((((((x4[t]-x1[t]))^2)+(((y4[t]-y1[t]))^2)+(((z4[t]-z1[t]))^2)))^3)])));


Fz1:=(((G (m1/(γ1z)) (m2/(γ2z)) (z2[t]-z1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((G (m1/(γ1z)) (m3/(γ3z)) (z3[t]-z1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((G  (m1/(γ1z)) (m4/(γ4z)) (z4[t]-z1[t]))/(Sqrt[((((((x4[t]-x1[t]))^2)+(((y4[t]-y1[t]))^2)+(((z4[t]-z1[t]))^2)))^3)]))+

((-q1*q2/(4π ε0 )(z2[t]-z1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((-q1*q3/(4π ε0 )(z3[t]-z1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((-q1*q4/(4π ε0 )(z4[t]-z1[t]))/(Sqrt[((((((x4[t]-x1[t]))^2)+(((y4[t]-y1[t]))^2)+(((z4[t]-z1[t]))^2)))^3)])));


Fx2:=(((G (m2/(γ2x)) (m1/(γ1x))  (x1[t]-x2[t]))/(Sqrt[((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)+(((z1[t]-z2[t]))^2)))^3)]))+((G  (m2/(γ2x)) (m3/(γ3x)) (x3[t]-x2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((G  (m2/(γ2x)) (m4/(γ4x)) (x4[t]-x2[t]))/(Sqrt[((((((x4[t]-x2[t]))^2)+(((y4[t]-y2[t]))^2)+(((z4[t]-z2[t]))^2)))^3)]))+

((-q2*q1/(4π ε0 ) (x1[t]-x2[t]))/(Sqrt[((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)+(((z1[t]-z2[t]))^2)))^3)]))+((-q2*q3/(4π ε0 )(x3[t]-x2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((-q2*q4/(4π ε0 )(x4[t]-x2[t]))/(Sqrt[((((((x4[t]-x2[t]))^2)+(((y4[t]-y2[t]))^2)+(((z4[t]-z2[t]))^2)))^3)])));


Fy2:=(((G (m2/(γ2y)) (m1/(γ1y))  (y1[t]-y2[t]))/(Sqrt[((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)+(((z1[t]-z2[t]))^2)))^3)]))+((G  (m2/(γ2y)) (m3/(γ3y)) (y3[t]-y2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((G  (m2/(γ2y)) (m4/(γ4y)) (y4[t]-y2[t]))/(Sqrt[((((((x4[t]-x2[t]))^2)+(((y4[t]-y2[t]))^2)+(((z4[t]-z2[t]))^2)))^3)]))+

((-q2*q1/(4π ε0 )(y1[t]-y2[t]))/(Sqrt[((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)+(((z1[t]-z2[t]))^2)))^3)]))+((-q2*q3/(4π ε0 )(y3[t]-y2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((-q2*q4/(4π ε0 )(y4[t]-y2[t]))/(Sqrt[((((((x4[t]-x2[t]))^2)+(((y4[t]-y2[t]))^2)+(((z4[t]-z2[t]))^2)))^3)])));


Fz2:=(((G (m2/(γ2z)) (m1/(γ1z))  (z1[t]-z2[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((G (m2/(γ2z)) (m3/(γ3z)) (z3[t]-z2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((G  (m2/(γ2z)) (m4/(γ4z)) (z4[t]-z2[t]))/(Sqrt[((((((x4[t]-x2[t]))^2)+(((y4[t]-y2[t]))^2)+(((z4[t]-z2[t]))^2)))^3)]))+

((-q2*q1/(4π ε0 )(z1[t]-z2[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((-q2*q3/(4π ε0 )(z3[t]-z2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((-q2*q4/(4π ε0 )(z4[t]-z2[t]))/(Sqrt[((((((x4[t]-x2[t]))^2)+(((y4[t]-y2[t]))^2)+(((z4[t]-z2[t]))^2)))^3)])));


Fx3:=(((G (m3/(γ3x)) (m1/(γ1x))  (x1[t]-x3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((G (m3/(γ3x)) (m2/(γ2x))  (x2[t]-x3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((G  (m3/(γ3x)) (m4/(γ4x)) (x4[t]-x3[t]))/(Sqrt[((((((x4[t]-x3[t]))^2)+(((y4[t]-y3[t]))^2)+(((z4[t]-z3[t]))^2)))^3)]))+

((-q3*q1/(4π ε0 )(x1[t]-x3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((-q3*q2/(4π ε0 )(x2[t]-x3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((-q3*q4/(4π ε0 )(x4[t]-x3[t]))/(Sqrt[((((((x4[t]-x3[t]))^2)+(((y4[t]-y3[t]))^2)+(((z4[t]-z3[t]))^2)))^3)])));


Fy3:=(((G (m3/(γ3y)) (m1/(γ1y))  (y1[t]-y3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((G (m3/(γ3y)) (m2/(γ2y))  (y2[t]-y3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((G  (m3/(γ3y)) (m4/(γ4y)) (y4[t]-y3[t]))/(Sqrt[((((((x4[t]-x3[t]))^2)+(((y4[t]-y3[t]))^2)+(((z4[t]-z3[t]))^2)))^3)]))+

((-q3*q1/(4π ε0 ) (y1[t]-y3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((-q3*q2/(4π ε0 ) (y2[t]-y3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((-q3*q4/(4π ε0 )(y4[t]-y3[t]))/(Sqrt[((((((x4[t]-x3[t]))^2)+(((y4[t]-y3[t]))^2)+(((z4[t]-z3[t]))^2)))^3)])));


Fz3:=(((G (m3/(γ3z)) (m1/(γ1z))  (z1[t]-z3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((G (m3/(γ3z)) (m2/(γ2z))  (z2[t]-z3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((G  (m3/(γ3z)) (m4/(γ4z)) (z4[t]-z3[t]))/(Sqrt[((((((x4[t]-x3[t]))^2)+(((y4[t]-y3[t]))^2)+(((z4[t]-z3[t]))^2)))^3)]))+

((-q3*q1/(4π ε0 )  (z1[t]-z3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((-q3*q2/(4π ε0 )  (z2[t]-z3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((-q3*q4/(4π ε0 )(z4[t]-z3[t]))/(Sqrt[((((((x4[t]-x3[t]))^2)+(((y4[t]-y3[t]))^2)+(((z4[t]-z3[t]))^2)))^3)])));


Fx4:=(((G (m4/(γ4x)) (m1/(γ1x))  (x1[t]-x4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((G (m4/(γ4x)) (m2/(γ2x))  (x2[t]-x4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((G  (m4/(γ4x)) (m3/(γ3x)) (x3[t]-x4[t]))/(Sqrt[((((((x3[t]-x4[t]))^2)+(((y3[t]-y4[t]))^2)+(((z3[t]-z4[t]))^2)))^3)]))+

((-q4*q1/(4π ε0 )(x1[t]-x4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((-q4*q2/(4π ε0 ) (x2[t]-x4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((-q4*q3/(4π ε0 ) (x3[t]-x4[t]))/(Sqrt[((((((x3[t]-x4[t]))^2)+(((y3[t]-y4[t]))^2)+(((z3[t]-z4[t]))^2)))^3)])));


Fy4:=(((G ((m4 )/(γ4y))(m1/(γ1y))  (y1[t]-y4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((G (m4/(γ4y)) (m2/(γ2y))  (y2[t]-y4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((G  (m4/(γ4y)) (m3/(γ3y)) (y3[t]-y4[t]))/(Sqrt[((((((x3[t]-x4[t]))^2)+(((y3[t]-y4[t]))^2)+(((z3[t]-z4[t]))^2)))^3)]))+

((-q4*q1/(4π ε0 ) (y1[t]-y4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((-q4*q2/(4π ε0 ) (y2[t]-y4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((-q4*q3/(4π ε0 )(y3[t]-y4[t]))/(Sqrt[((((((x3[t]-x4[t]))^2)+(((y3[t]-y4[t]))^2)+(((z3[t]-z4[t]))^2)))^3)])));


Fz4:=(((G (m4/(γ4z)) (m1/(γ1z))  (z1[t]-z4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((G (m4/(γ4z)) (m2/(γ2z))  (z2[t]-z4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((G  (m4/(γ4z)) (m3/(γ3z)) (z3[t]-z4[t]))/(Sqrt[((((((x3[t]-x4[t]))^2)+(((y3[t]-y4[t]))^2)+(((z3[t]-z4[t]))^2)))^3)]))+

((-q4*q1/(4π ε0 )(z1[t]-z4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((-q4*q2/(4π ε0 )(z2[t]-z4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((-q4*q3/(4π ε0 )(z3[t]-z4[t]))/(Sqrt[((((((x3[t]-x4[t]))^2)+(((y3[t]-y4[t]))^2)+(((z3[t]-z4[t]))^2)))^3)])));


m1=M1; m2=M2; m3=M3; m4=M4;q1=Q1; q2= Q2; q3=Q3; q4=Q4;
x1x=P1x; x2x=P2x; x3x=P3x; x4x=P4x;y1y=P1y; y2y=P2y; y3y=P3y; y4y=P4y;z1z=P1z;z2z=P2z; z3z=P3z; z4z=P4z;
vx1x=v1x; vx2x=v2x; vx3x=v3x; vx4x=v4x;vy1y=v1y; vy2y=v2y; vy3y=v3y; vy4y=v4y;vz1z=v1z;vz2z=v2z; vz3z=v3z; vz4z=v4z;


(* Container *)
Funktion[{{x1x_,y1y_,z1z_},{x2x_,y2y_,z2z_},{x3x_,y3y_,z3z_},{x4x_,y4y_,z4z_}},
{{vx1x_,vy1y_,vz1z_},{vx2x_,vy2y_,vz2z_},{vx3x_,vy3y_,vz3z_},{vx4x_,vy4y_,vz4z_}},
{m1_,m2_,m3_,m4_},
{q1_,q2_,q3_,q4_},
T_,plotType:("x"|"v"),plotOptions___]:=

Module[{nds,Tmax,funcToPlot},

(* DGL *)
nds=NDSolve[{

x1'[t]==vx1[t],y1'[t]==vy1[t],z1'[t]==vz1[t],x2'[t]==vx2[t],y2'[t]==vy2[t],z2'[t]==vz2[t],
x3'[t]==vx3[t],y3'[t]==vy3[t],z3'[t]==vz3[t],x4'[t]==vx4[t],y4'[t]==vy4[t],z4'[t]==vz4[t],

(* Beschleunigung *)

vx1'[t]==((Fx1-((vx1[t](Fx1 vx1[t]+Fy1 vy1[t]+Fz1 vz1[t]))/(c^2)))/(Sqrt[((((m1(Sqrt[(((vx1[t]))^2)+(((vy1[t]))^2)+(((vz1[t]))^2)])))^2)/((c^2)(1-(((((vx1[t]))^2)+(((vy1[t]))^2)+(((vz1[t]))^2))/(c^2)))))+(m1^2)])), vy1'[t]==((Fy1-((vy1[t](Fx1 vx1[t]+Fy1 vy1[t]+Fz1 vz1[t]))/(c^2)))/(Sqrt[((((m1(Sqrt[(((vx1[t]))^2)+(((vy1[t]))^2)+(((vz1[t]))^2)])))^2)/((c^2)(1-(((((vx1[t]))^2)+(((vy1[t]))^2)+(((vz1[t]))^2))/(c^2)))))+(m1^2)])), vz1'[t]==((Fz1-((vz1[t](Fx1 vx1[t]+Fy1 vy1[t]+Fz1 vz1[t]))/(c^2)))/(Sqrt[((((m1(Sqrt[(((vx1[t]))^2)+(((vy1[t]))^2)+(((vz1[t]))^2)])))^2)/((c^2)(1-(((((vx1[t]))^2)+(((vy1[t]))^2)+(((vz1[t]))^2))/(c^2)))))+(m1^2)])),

vx2'[t]==((Fx2-((vx2[t](Fx2 vx2[t]+Fy2 vy2[t]+Fz2 vz2[t]))/(c^2)))/(Sqrt[((((m2(Sqrt[(((vx2[t]))^2)+(((vy2[t]))^2)+(((vz2[t]))^2)])))^2)/((c^2)(1-(((((vx2[t]))^2)+(((vy2[t]))^2)+(((vz2[t]))^2))/(c^2)))))+(m2^2)])), vy2'[t]==((Fy2-((vy2[t](Fx2 vx2[t]+Fy2 vy2[t]+Fz2 vz2[t]))/(c^2)))/(Sqrt[((((m2(Sqrt[(((vx2[t]))^2)+(((vy2[t]))^2)+(((vz2[t]))^2)])))^2)/((c^2)(1-(((((vx2[t]))^2)+(((vy2[t]))^2)+(((vz2[t]))^2))/(c^2)))))+(m2^2)])), vz2'[t]==((Fz2-((vz2[t](Fx2 vx2[t]+Fy2 vy2[t]+Fz2 vz2[t]))/(c^2)))/(Sqrt[((((m2(Sqrt[(((vx2[t]))^2)+(((vy2[t]))^2)+(((vz2[t]))^2)])))^2)/((c^2)(1-(((((vx2[t]))^2)+(((vy2[t]))^2)+(((vz2[t]))^2))/(c^2)))))+(m2^2)])),

vx3'[t]==((Fx3-((vx3[t](Fx3 vx3[t]+Fy3 vy3[t]+Fz3 vz3[t]))/(c^2)))/(Sqrt[((((m3(Sqrt[(((vx3[t]))^2)+(((vy3[t]))^2)+(((vz3[t]))^2)])))^2)/((c^2)(1-(((((vx3[t]))^2)+(((vy3[t]))^2)+(((vz3[t]))^2))/(c^2)))))+(m3^2)])), vy3'[t]==((Fy3-((vy3[t](Fx3 vx3[t]+Fy3 vy3[t]+Fz3 vz3[t]))/(c^2)))/(Sqrt[((((m3(Sqrt[(((vx3[t]))^2)+(((vy3[t]))^2)+(((vz3[t]))^2)])))^2)/((c^2)(1-(((((vx3[t]))^2)+(((vy3[t]))^2)+(((vz3[t]))^2))/(c^2)))))+(m3^2)])), vz3'[t]==((Fz3-((vz3[t](Fx3 vx3[t]+Fy3 vy3[t]+Fz3 vz3[t]))/(c^2)))/(Sqrt[((((m3(Sqrt[(((vx3[t]))^2)+(((vy3[t]))^2)+(((vz3[t]))^2)])))^2)/((c^2)(1-(((((vx3[t]))^2)+(((vy3[t]))^2)+(((vz3[t]))^2))/(c^2)))))+(m3^2)])),

vx4'[t]==((Fx4-((vx4[t](Fx4 vx4[t]+Fy4 vy4[t]+Fz4 vz4[t]))/(c^2)))/(Sqrt[((((m4(Sqrt[(((vx4[t]))^2)+(((vy4[t]))^2)+(((vz4[t]))^2)])))^2)/((c^2)(1-(((((vx4[t]))^2)+(((vy4[t]))^2)+(((vz4[t]))^2))/(c^2)))))+(m4^2)])), vy4'[t]==((Fy4-((vy4[t](Fx4 vx4[t]+Fy4 vy4[t]+Fz4 vz4[t]))/(c^2)))/(Sqrt[((((m4(Sqrt[(((vx4[t]))^2)+(((vy4[t]))^2)+(((vz4[t]))^2)])))^2)/((c^2)(1-(((((vx4[t]))^2)+(((vy4[t]))^2)+(((vz4[t]))^2))/(c^2)))))+(m4^2)])), vz4'[t]==((Fz4-((vz4[t](Fx4 vx4[t]+Fy4 vy4[t]+Fz4 vz4[t]))/(c^2)))/(Sqrt[((((m4(Sqrt[(((vx4[t]))^2)+(((vy4[t]))^2)+(((vz4[t]))^2)])))^2)/((c^2)(1-(((((vx4[t]))^2)+(((vy4[t]))^2)+(((vz4[t]))^2))/(c^2)))))+(m4^2)])),

x1[0]==x1x,y1[0]==y1y,z1[0]==z1z,x2[0]==x2x,y2[0]==y2y,z2[0]==z2z,
x3[0]==x3x,y3[0]==y3y,z3[0]==z3z,x4[0]==x4x,y4[0]==y4y,z4[0]==z4z,

vx1[0]==vx1x,vy1[0]==vy1y,vz1[0]==vz1z,vx2[0]==vx2x,vy2[0]==vy2y,vz2[0]==vz2z,
vx3[0]==vx3x,vy3[0]==vy3y,vz3[0]==vz3z,vx4[0]==vx4x,vy4[0]==vy4y,vz4[0]==vz4z},

{x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
vx1,vx2,vx3,vx4,vy1,vy2,vy3,vy4,vz1,vz2,vz3,vz4},

{t,0,T},

MaxSteps->Infinity,
Method->{"ImplicitRungeKutta",
"DifferenceOrder"->20},
InterpolationOrder->All];

(* Solve *)
If[Head[nds]=!=NDSolve,Tmax=nds[[1,1,2,1,1,2]];
funcToPlot=If[plotType==="x",
{{x1[t],y1[t],z1[t]},{x2[t],y2[t],z2[t]},{x3[t],y3[t],z3[t]},{x4[t],y4[t],z4[t]}},
{{vx1[t],vy1[t],vz1[t]},{vx2[t],vy2[t],vz2[t]},{vx3[t],vy3[t],vz3[t]},{vx4[t],vy4[t],vz4[t]}}]/.nds[[1]];

(* Plot *)
ParametricPlot3D[Evaluate[funcToPlot],

(* Schweiflänge *)
{t,Tmax-Tmin,Tmax},

(* Farben *)
PlotStyle->{{Red,Thick},{Blue,Thick},{Green,Thick},{Orange,Thick}},

(* Plot Range *)
PlotRange->{{-2 Au,2 Au},{-2 Au,2 Au},{-2 Au,2 Au}},
AspectRatio->1,Ticks->None,
MaxRecursion->ControlActive[3,Infinity],
Axes->True,plotOptions],Text["Yukterez Mod."]]]//Quiet

(* Manipulate *)
plot[a_, b_, c_]:=Quiet[Manipulate[Show[
Funktion[{
{P1x,P1y,P1z},
{P2x,P2y,P2z},
{P3x,P3y,P3z},
{P4x,P4y,P4z}},
{{v1x,v1y,v1z},
{v2x,v2y,v2z},
{v3x,v3y,v3z},
{v4x,v4y,v4z}},
{M1,M2,M3,M4},
{Q1,Q2,Q3,Q4},
T,xv,ImageSize->{380, 380}],

(* Startpositionen *)
Graphics3D[{
{LightRed,Point[{P1x,P1y,P1z}]},
{LightBlue,Point[{P2x,P2y,P2z}]},
{LightGreen,Point[{P3x,P3y,P3z}]},
{LightOrange,Point[{P4x,P4y,P4z}]}}],

(* Perspektive *)
ViewPoint->{a,b,c}],
{{xv,"x","Modus"},{"x"->"Position"}},

(* Zeitregler *)
{{T,Tmin,"Time"},Tmin,τ},ControlPlacement->{Bottom}]];

p1=plot[0,0,4]
p2=plot[0,4,0]
p3=plot[4,4,4]

Quit[]