Für den vollständigen Artikel geht es hier entlang; für die Diskussion klick hier. n-Body Solver: (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* | Mathematica Syntax | yukterez.net | n-Body calculator for disks with height | Version 1 | *) (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) G  = 1;                                                               (* Gravitationskonstante *) d  = 0.002;                                         (* Abstand der Scheibenpartikel zueinander *) n  = Floor[2π/d];                                                        (* Abstandsgleichheit *) я1 = d;                                                                 (* Scheibeninnenradius *) я2 = 1;                                                                 (* Scheibenaußenradius *) h  = 0.022;                                                                    (* Scheibenhöhe *) fx[x_, y_, z_] :=                                                                                 Flatten[Table[Table[Table[(G M (r Sin[φ]-x))/Sqrt[((r Sin[φ]-x)^2+(r Cos[φ]-y)^2+(Z-z)^2)^3],     {φ, 0, 2π-2π/r/n, 2π/r/n}], {r, я1, я2, d}], {Z, -h/2, h/2, d}]]                                  gx[x_, y_, z_] := Total[-fx[x, y, z]]                          (* Fallbeschleunigung entlang x *) fy[x_, y_, z_] :=                                                                                 Flatten[Table[Table[Table[(G M (r Cos[φ]-y))/Sqrt[((r Sin[φ]-x)^2+(r Cos[φ]-y)^2+(Z-z)^2)^3],     {φ, 0, 2π-2π/r/n, 2π/r/n}], {r, я1, я2, d}], {Z, -h/2, h/2, d}]]                                  gy[x_, y_, z_] := Total[-fy[x, y, z]]                          (* Fallbeschleunigung entlang y *) fz[x_, y_, z_] :=                                                                                 Flatten[Table[Table[Table[(G M (Z-z))/Sqrt[((r Sin[φ]-x)^2+(r Cos[φ]-y)^2+(Z-z)^2)^3],            {φ, 0, 2π-2π/r/n, 2π/r/n}], {r, я1, я2, d}], {Z, -h/2, h/2, d}]]                                  gz[x_, y_, z_] := Total[-fz[x, y, z]]                          (* Fallbeschleunigung entlang z *) l = Length[fx[0, 0, 0]]                   (* Anzahl der Partikel aus denen die Scheibe besteht *) M = 1/l;                                                       (* Masse der einzelnen Partikel *) gχ = Interpolation[ParallelTable[{x, gx[x, 0, 0]}, {x, 0.01, 2, 0.01}]]                           Plot[gχ[x], {x, 0, 2}, Frame->True, ImageSize->640, PlotRange->{{0, 2}, {0, 3.5}},                AspectRatio->1/3, ImagePadding->{{25, 1}, {15, 5}}, GridLines->{{1}, {1}}]