Here are the best pictures I have made. The C code uses fourth order Runge-Kutta with variable Δt. See Weisstein’s note on Runge-Kutta for what you need to know to write the code. (See at least these pictures for clues as to why it works.) Here are comments on the accuracy of the method. 259818 steps were used. This was computed with the Pentium’s 64 significant bits. This Java code (explained here) reads data files produced by the C code and produces .png files which are shown below.

Each picture is for 8 time units. White dots on the curve are where the time has an integral value. Labels for times 35 thru 40 are confusing. At t = 37 and t = 38 the 5 mass (yellow) lands on the track of the 3 mass (turquois). At t = 47 the 5 mass is just back to where it was at about t = 40.6. The greatest numerical distress occurs about 15.7 which event is notable in the picture. The 3 mass is off screen between t = 49 and t = 57 after which it makes its last encounter with the other two which have already begun their final perpetual pas de deux. It passes between the other two in such a phase so as to extract enough of their potential energy to be ejected from the system.

Legend:
masscolor
3turquois
4magenta
5yellow
0<=t<8 at t = 1.8794
ke = 2035
Δt = 4*10−7

at t = 3.8
ke = 324
Δt = 6.2*10−6

at t = 6.898
ke = 140
Δt = 1.3*10−5

8<=t<16 at t = 8.76
ke = 2341
Δt = 3.2*10−7

at t = 15.82992
ke = 48249 (Wow!)
Δt = 3.4*10−9

16<=t<24 at t = 22.965
ke = 1137
Δt = 9.3*10−7
24<=t<32 at t = 29.8015
ke = 7153
Δt = 5.9*10−8
32<=t<40 at t = 33.6715
ke = 213
Δt = 7.6*10−6
40<=t<48 at t = 41.209
ke = 619
Δt = 2.4*10−6

at t = 46.015
ke = 596
Δt = 2.5*10−6

at t = 47.371
ke = 1189
Δt = 8.7*10−7

48<=t<56 at t = 52.116
ke = 1513
Δt = 6*10−7

etc.

56<=t<62 at t = 59.777
ke = 2662
Δt = 6.2*10−8

at t = 60.6335
ke = 3067
Δt = 6.2*10−8

at t = 61.494
ke = 3126
Δt = 3.1*10−8

Accuracy