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}}]