@Carl Love  Thanks a lot for your code but I always have this message:

Error, (in plots/animate) input is not a valid dsolve/numeric solution.

I don't know why?!

restart:
with(linalg):
with(DEtools):

#Positions des éléments de la galaxie

Position1 := [x1(t), y1(t), z1(t)]:
Position2 := [x2(t), y2(t), z2(t)]:
Position3 := [x3(t), y3(t), z3(t)]:
Position4 := [x4(t), y4(t), z4(t)]:
Position5 := [x5(t), y5(t), z5(t)]:
Position6 := [x6(t), y6(t), z6(t)]:
Position7 := [x7(t), y7(t), z7(t)]:
Position8 := [x8(t), y8(t), z8(t)]:
Position9 := [x9(t), y9(t), z9(t)]:
PositionA := [xA(t), yA(t), zA(t)]:

#Distance entre les étoiles (unité de distance D)
r32 := sqrt((x3(t) - x2(t))^2 + (y3(t) - y2(t))^2 + (z3(t) - z2(t))^2 + 0.001):
r31 := sqrt((x3(t) - x1(t))^2 + (y3(t) - y1(t))^2 + (z3(t) - z1(t))^2 + 0.001):
r12 := sqrt((x1(t) - x2(t))^2 + (y1(t) - y2(t))^2 + (z1(t) - z2(t))^2 + 0.001):
r42 := sqrt((x4(t) - x2(t))^2 + (y4(t) - y2(t))^2 + (z4(t) - z2(t))^2 + 0.001):
r41 := sqrt((x4(t) - x1(t))^2 + (y4(t) - y1(t))^2 + (z4(t) - z1(t))^2 + 0.001):
r43 := sqrt((x4(t) - x3(t))^2 + (y4(t) - y3(t))^2 + (z4(t) - z3(t))^2 + 0.001):
r51 := sqrt((x5(t) - x1(t))^2 + (y5(t) - y1(t))^2 + (z5(t) - z1(t))^2 + 0.001):
r52 := sqrt((x5(t) - x2(t))^2 + (y5(t) - y2(t))^2 + (z5(t) - z2(t))^2 + 0.001):
r53 := sqrt((x5(t) - x3(t))^2 + (y5(t) - y3(t))^2 + (z5(t) - z3(t))^2 + 0.001):
r54 := sqrt((x5(t) - x4(t))^2 + (y5(t) - y4(t))^2 + (z5(t) - z4(t))^2 + 0.001):
r61 := sqrt((x6(t) - x1(t))^2 + (y6(t) - y1(t))^2 + (z6(t) - z1(t))^2 + 0.001):
r62 := sqrt((x6(t) - x2(t))^2 + (y6(t) - y2(t))^2 + (z6(t) - z2(t))^2 + 0.001):
r63 := sqrt((x6(t) - x3(t))^2 + (y6(t) - y3(t))^2 + (z6(t) - z3(t))^2 + 0.001):
r64 := sqrt((x6(t) - x4(t))^2 + (y6(t) - y4(t))^2 + (z6(t) - z4(t))^2 + 0.001):
r65 := sqrt((x6(t) - x5(t))^2 + (y6(t) - y5(t))^2 + (z6(t) - z5(t))^2 + 0.001):
r71 := sqrt((x7(t) - x1(t))^2 + (y7(t) - y1(t))^2 + (z7(t) - z1(t))^2 + 0.001):
r72 := sqrt((x7(t) - x2(t))^2 + (y7(t) - y2(t))^2 + (z7(t) - z2(t))^2 + 0.001):
r73 := sqrt((x7(t) - x3(t))^2 + (y7(t) - y3(t))^2 + (z7(t) - z3(t))^2 + 0.001):
r74 := sqrt((x7(t) - x4(t))^2 + (y7(t) - y4(t))^2 + (z7(t) - z4(t))^2 + 0.001):
r75 := sqrt((x7(t) - x5(t))^2 + (y7(t) - y5(t))^2 + (z7(t) - z5(t))^2 + 0.001):
r76 := sqrt((x7(t) - x6(t))^2 + (y7(t) - y6(t))^2 + (z7(t) - z6(t))^2 + 0.001):
r81 := sqrt((x8(t) - x1(t))^2 + (y8(t) - y1(t))^2 + (z8(t) - z1(t))^2 + 0.001):
r82 := sqrt((x8(t) - x2(t))^2 + (y8(t) - y2(t))^2 + (z8(t) - z2(t))^2 + 0.001):
r83 := sqrt((x8(t) - x3(t))^2 + (y8(t) - y3(t))^2 + (z8(t) - z3(t))^2 + 0.001):
r84 := sqrt((x8(t) - x4(t))^2 + (y8(t) - y4(t))^2 + (z8(t) - z4(t))^2 + 0.001):
r85 := sqrt((x8(t) - x5(t))^2 + (y8(t) - y5(t))^2 + (z8(t) - z5(t))^2 + 0.001):
r86 := sqrt((x8(t) - x6(t))^2 + (y8(t) - y6(t))^2 + (z8(t) - z6(t))^2 + 0.001):
r87 := sqrt((x8(t) - x7(t))^2 + (y8(t) - y7(t))^2 + (z8(t) - z7(t))^2 + 0.001):
r91 := sqrt((x9(t) - x1(t))^2 + (y9(t) - y1(t))^2 + (z9(t) - z1(t))^2 + 0.001):
r92 := sqrt((x9(t) - x2(t))^2 + (y9(t) - y2(t))^2 + (z9(t) - z2(t))^2 + 0.001):
r93 := sqrt((x9(t) - x3(t))^2 + (y9(t) - y3(t))^2 + (z9(t) - z3(t))^2 + 0.001):
r94 := sqrt((x9(t) - x4(t))^2 + (y9(t) - y4(t))^2 + (z9(t) - z4(t))^2 + 0.001):
r95 := sqrt((x9(t) - x5(t))^2 + (y9(t) - y5(t))^2 + (z9(t) - z5(t))^2 + 0.001):
r96 := sqrt((x9(t) - x6(t))^2 + (y9(t) - y6(t))^2 + (z9(t) - z6(t))^2 + 0.001):
r97 := sqrt((x9(t) - x7(t))^2 + (y9(t) - y7(t))^2 + (z9(t) - z7(t))^2 + 0.001):
r98 := sqrt((x9(t) - x8(t))^2 + (y9(t) - y8(t))^2 + (z9(t) - z8(t))^2 + 0.001):
rA1 := sqrt((xA(t) - x1(t))^2 + (yA(t) - y1(t))^2 + (zA(t) - z1(t))^2 + 0.001):
rA2 := sqrt((xA(t) - x2(t))^2 + (yA(t) - y2(t))^2 + (zA(t) - z2(t))^2 + 0.001):
rA3 := sqrt((xA(t) - x3(t))^2 + (yA(t) - y3(t))^2 + (zA(t) - z3(t))^2 + 0.001):
rA4 := sqrt((xA(t) - x4(t))^2 + (yA(t) - y4(t))^2 + (zA(t) - z4(t))^2 + 0.001):
rA5 := sqrt((xA(t) - x5(t))^2 + (yA(t) - y5(t))^2 + (zA(t) - z5(t))^2 + 0.001):
rA6 := sqrt((xA(t) - x6(t))^2 + (yA(t) - y6(t))^2 + (zA(t) - z6(t))^2 + 0.001):
rA7 := sqrt((xA(t) - x7(t))^2 + (yA(t) - y7(t))^2 + (zA(t) - z7(t))^2 + 0.001):
rA8 := sqrt((xA(t) - x8(t))^2 + (yA(t) - y8(t))^2 + (zA(t) - z8(t))^2 + 0.001):
rA9 := sqrt((xA(t) - x9(t))^2 + (yA(t) - y9(t))^2 + (zA(t) - z9(t))^2 + 0.001):

#Vitesse des éléments de la galaxie (en D/T)
Vitesse1 := [(D(x1))(t), (D(y1))(t), (D(z1))(t)]:
Vitesse2 := [(D(x2))(t), (D(y2))(t), (D(z2))(t)]:
Vitesse3 := [(D(x3))(t), (D(y3))(t), (D(z3))(t)]:
Vitesse4 := [(D(x4))(t), (D(y4))(t), (D(z4))(t)]:
Vitesse5 := [(D(x5))(t), (D(y5))(t), (D(z5))(t)]:
Vitesse6 := [(D(x6))(t), (D(y6))(t), (D(z6))(t)]:
Vitesse7 := [(D(x7))(t), (D(y7))(t), (D(z7))(t)]:
Vitesse8 := [(D(x8))(t), (D(y8))(t), (D(z8))(t)]:
Vitesse9 := [(D(x9))(t), (D(y9))(t), (D(z9))(t)]:
VitesseA := [(D(xA))(t), (D(yA))(t), (D(zA))(t)]:

#Masse des étoiles de la galaxie (unité de masse M)
Masse1 := 100.0:
Masse2 := 0.00001667:
Masse3 := 0.0002459:
Masse4 := 0.0003015:
Masse5 := 0.00003237:
Masse6 := 1.0:
Masse7 := 0.028667:
Masse8 := 0.000437:
Masse9 := 0.00520202:
MasseA := 5.0:

#Force gravitationnelle (unité de force N)
Fg32 := -Masse2*(Position3 - Position2)/((r32)^3):
Fg23 := -Masse3*(Position2 - Position3)/((r32)^3):
Fg13 := -Masse3*(Position1 - Position3)/((r31)^3):
Fg31 := -Masse1*(Position3 - Position1)/((r31)^3):
Fg21 := -Masse1*(Position2 - Position1)/((r12)^3):
Fg12 := -Masse2*(Position1 - Position2)/((r12)^3):
Fg41 := -Masse1*(Position4 - Position1)/((r41)^3):
Fg14 := -Masse4*(Position1 - Position4)/((r41)^3):
Fg42 := -Masse2*(Position4 - Position2)/((r42)^3):
Fg24 := -Masse4*(Position2 - Position4)/((r42)^3):
Fg43 := -Masse3*(Position4 - Position3)/((r43)^3):
Fg34 := -Masse4*(Position3 - Position4)/((r43)^3):
Fg51 := -Masse1*(Position5 - Position1)/((r51)^3):
Fg15 := -Masse5*(Position1 - Position5)/((r51)^3):
Fg52 := -Masse2*(Position5 - Position2)/((r52)^3):
Fg25 := -Masse5*(Position2 - Position5)/((r52)^3):
Fg53 := -Masse3*(Position5 - Position3)/((r53)^3):
Fg35 := -Masse5*(Position3 - Position5)/((r53)^3):
Fg54 := -Masse4*(Position5 - Position4)/((r54)^3):
Fg45 := -Masse5*(Position4 - Position5)/((r54)^3):
Fg61 := -Masse1*(Position6 - Position1)/((r61)^3):
Fg16 := -Masse6*(Position1 - Position6)/((r61)^3):
Fg62 := -Masse2*(Position6 - Position2)/((r62)^3):
Fg26 := -Masse6*(Position2 - Position6)/((r62)^3):
Fg63 := -Masse3*(Position6 - Position3)/((r63)^3):
Fg36 := -Masse6*(Position3 - Position6)/((r63)^3):
Fg64 := -Masse4*(Position6 - Position4)/((r64)^3):
Fg46 := -Masse6*(Position4 - Position6)/((r64)^3):
Fg65 := -Masse5*(Position6 - Position5)/((r65)^3):
Fg56 := -Masse6*(Position5 - Position6)/((r65)^3):
Fg71 := -Masse1*(Position7 - Position1)/((r71)^3):
Fg17 := -Masse7*(Position1 - Position7)/((r71)^3):
Fg72 := -Masse2*(Position7 - Position2)/((r72)^3):
Fg27 := -Masse7*(Position2 - Position7)/((r72)^3):
Fg73 := -Masse3*(Position7 - Position3)/((r73)^3):
Fg37 := -Masse7*(Position3 - Position7)/((r73)^3):
Fg74 := -Masse4*(Position7 - Position4)/((r74)^3):
Fg47 := -Masse7*(Position4 - Position7)/((r74)^3):
Fg75 := -Masse5*(Position7 - Position5)/((r75)^3):
Fg57 := -Masse7*(Position5 - Position7)/((r75)^3):
Fg76 := -Masse6*(Position7 - Position6)/((r76)^3):
Fg67 := -Masse7*(Position6 - Position7)/((r76)^3):
Fg81 := -Masse1*(Position8 - Position1)/((r81)^3):
Fg18 := -Masse8*(Position1 - Position8)/((r81)^3):
Fg82 := -Masse2*(Position8 - Position2)/((r82)^3):
Fg28 := -Masse8*(Position2 - Position8)/((r82)^3):
Fg83 := -Masse3*(Position8 - Position3)/((r83)^3):
Fg38 := -Masse8*(Position3 - Position8)/((r83)^3):
Fg84 := -Masse4*(Position8 - Position4)/((r84)^3):
Fg48 := -Masse8*(Position4 - Position8)/((r84)^3):
Fg85 := -Masse5*(Position8 - Position5)/((r85)^3):
Fg58 := -Masse8*(Position5 - Position8)/((r85)^3):
Fg86 := -Masse6*(Position8 - Position6)/((r86)^3):
Fg68 := -Masse8*(Position6 - Position8)/((r86)^3):
Fg87 := -Masse7*(Position8 - Position7)/((r87)^3):
Fg78 := -Masse8*(Position7 - Position8)/((r87)^3):
Fg91 := -Masse1*(Position9 - Position1)/((r91)^3):
Fg19 := -Masse9*(Position1 - Position9)/((r91)^3):
Fg92 := -Masse2*(Position9 - Position2)/((r92)^3):
Fg29 := -Masse9*(Position2 - Position9)/((r92)^3):
Fg93 := -Masse3*(Position9 - Position3)/((r93)^3):
Fg39 := -Masse9*(Position3 - Position9)/((r93)^3):
Fg94 := -Masse4*(Position9 - Position4)/((r94)^3):
Fg49 := -Masse9*(Position4 - Position9)/((r94)^3):
Fg95 := -Masse5*(Position9 - Position5)/((r95)^3):
Fg59 := -Masse9*(Position5 - Position9)/((r95)^3):
Fg96 := -Masse6*(Position9 - Position6)/((r96)^3):
Fg69 := -Masse9*(Position6 - Position9)/((r96)^3):
Fg97 := -Masse7*(Position9 - Position7)/((r97)^3):
Fg79 := -Masse9*(Position7 - Position9)/((r97)^3):
Fg98 := -Masse8*(Position9 - Position8)/((r98)^3):
Fg89 := -Masse9*(Position8 - Position9)/((r98)^3):
FgA1 := -Masse1*(PositionA - Position1)/((rA1)^3):
Fg1A := -MasseA*(Position1 - PositionA)/((rA1)^3):
FgA2 := -Masse2*(PositionA - Position2)/((rA2)^3):
Fg2A := -MasseA*(Position2 - PositionA)/((rA2)^3):
FgA3 := -Masse3*(PositionA - Position3)/((rA3)^3):
Fg3A := -MasseA*(Position3 - PositionA)/((rA3)^3):
FgA4 := -Masse4*(PositionA - Position4)/((rA4)^3):
Fg4A := -MasseA*(Position4 - PositionA)/((rA4)^3):
FgA5 := -Masse5*(PositionA - Position5)/((rA5)^3):
Fg5A := -MasseA*(Position5 - PositionA)/((rA5)^3):
FgA6 := -Masse6*(PositionA - Position6)/((rA6)^3):
Fg6A := -MasseA*(Position6 - PositionA)/((rA6)^3):
FgA7 := -Masse7*(PositionA - Position7)/((rA7)^3):
Fg7A := -MasseA*(Position7 - PositionA)/((rA7)^3):
FgA8 := -Masse8*(PositionA - Position8)/((rA8)^3):
Fg8A := -MasseA*(Position8 - PositionA)/((rA8)^3):
FgA9 := -Masse9*(PositionA - Position9)/((rA9)^3):
Fg9A := -MasseA*(Position9 - PositionA)/((rA9)^3):

#Force exercée sur chaque étoile (N)
Force1:= Fg12+Fg13+Fg14+Fg15+Fg16+Fg17+Fg18+Fg19+Fg1A:
Force2:= Fg21+Fg23+Fg24+Fg25+Fg26+Fg27+Fg28+Fg29+Fg2A:
Force3:= Fg31+Fg32+Fg34+Fg35+Fg36+Fg37+Fg38+Fg39+Fg3A:
Force4:= Fg41+Fg42+Fg43+Fg45+Fg46+Fg47+Fg48+Fg49+Fg4A:
Force5:= Fg51+Fg52+Fg53+Fg54+Fg56+Fg57+Fg58+Fg59+Fg5A:
Force6:= Fg61+Fg62+Fg63+Fg64+Fg65+Fg67+Fg68+Fg69+Fg6A:
Force7:= Fg71+Fg72+Fg73+Fg74+Fg75+Fg76+Fg78+Fg79+Fg7A:
Force8:= Fg81+Fg82+Fg83+Fg84+Fg85+Fg86+Fg87+Fg89+Fg8A:
Force9:= Fg91+Fg92+Fg93+Fg94+Fg95+Fg96+Fg97+Fg98+Fg9A:
ForceA:= FgA1+FgA2+FgA3+FgA4+FgA5+FgA6+FgA7+FgA8+FgA9:

#Produits scalaire...(???)
F1x := innerprod([1, 0, 0], Force1):
F1y := innerprod([0, 1, 0], Force1):
F1z := innerprod([0, 0, 1], Force1):
F2x := innerprod([1, 0, 0], Force2):
F2y := innerprod([0, 1, 0], Force2):
F2z := innerprod([0, 0, 1], Force2):
F3x := innerprod([1, 0, 0], Force3):
F3y := innerprod([0, 1, 0], Force3):
F3z := innerprod([0, 0, 1], Force3):
F4x := innerprod([1, 0, 0], Force4):
F4y := innerprod([0, 1, 0], Force4):
F4z := innerprod([0, 0, 1], Force4):
F5x := innerprod([1, 0, 0], Force5):
F5y := innerprod([0, 1, 0], Force5):
F5z := innerprod([0, 0, 1], Force5):
F6x := innerprod([1, 0, 0], Force6):
F6y := innerprod([0, 1, 0], Force6):
F6z := innerprod([0, 0, 1], Force6):
F7x := innerprod([1, 0, 0], Force7):
F7y := innerprod([0, 1, 0], Force7):
F7z := innerprod([0, 0, 1], Force7):
F8x := innerprod([1, 0, 0], Force8):
F8y := innerprod([0, 1, 0], Force8):
F8z := innerprod([0, 0, 1], Force8):
F9x := innerprod([1, 0, 0], Force9):
F9y := innerprod([0, 1, 0], Force9):
F9z := innerprod([0, 0, 1], Force9):
FAx := innerprod([1, 0, 0], ForceA):
FAy := innerprod([0, 1, 0], ForceA):
FAz := innerprod([0, 0, 1], ForceA):

#Temps de la simulation (unité de temps T)
TempsInitial := -100:
TempsFinal := 1600:

#Équations de mes composantes
eq1x := (D(D(x1)))(t) = F1x:
eq1y := (D(D(y1)))(t) = F1y:
eq1z := (D(D(z1)))(t) = F1z:
eq2x := (D(D(x2)))(t) = F2x:
eq2y := (D(D(y2)))(t) = F2y:
eq2z := (D(D(z2)))(t) = F2z:
eq3x := (D(D(x3)))(t) = F3x:
eq3y := (D(D(y3)))(t) = F3y:
eq3z := (D(D(z3)))(t) = F3z:
eq4x := (D(D(x4)))(t) = F4x:
eq4y := (D(D(y4)))(t) = F4y:
eq4z := (D(D(z4)))(t) = F4z:
eq5x := (D(D(x5)))(t) = F5x:
eq5y := (D(D(y5)))(t) = F5y:
eq5z := (D(D(z5)))(t) = F5z:
eq6x := (D(D(x6)))(t) = F6x:
eq6y := (D(D(y6)))(t) = F6y:
eq6z := (D(D(z6)))(t) = F6z:
eq7x := (D(D(x7)))(t) = F7x:
eq7y := (D(D(y7)))(t) = F7y:
eq7z := (D(D(z7)))(t) = F7z:
eq8x := (D(D(x8)))(t) = F8x:
eq8y := (D(D(y8)))(t) = F8y:
eq8z := (D(D(z8)))(t) = F8z:
eq9x := (D(D(x9)))(t) = F9x:
eq9y := (D(D(y9)))(t) = F9y:
eq9z := (D(D(z9)))(t) = F9z:
eqAx := (D(D(xA)))(t) = FAx:
eqAy := (D(D(yA)))(t) = FAy:
eqAz := (D(D(zA)))(t) = FAz:

#Conditions initiales de la simulation

x10 := 0.0:
y10 := 0.0:
z10 := 0.0:
V1x0 := 0.0:
V1y0 := 0.02615894084:
V1z0 := 0.0:

x20 := 1.17:
y20 := 0.0:
z20 := 0.0:
V2x0 := 0.0:
V2y0 := -9.245:
V2z0 := 0.0:

x30 := 2.169:
y30 := 0.0:
z30 := 0.0:
V3x0 := 0.0:
V3y0 := -6.79:
V3z0 := 0.0:

x40 := 3.0:
y40 := 0.0:
z40 := 0.0:
V4x0 := 0.0:
V4y0 := -5.774:
V4z0 := 0.0:

x50 := 4.56:
y50 := 0.0:
z50 := 0.0:
V5x0 := 0.0:
V5y0 := -4.683:
V5z0 := 0.0:

x60 := 15.6:
y60 := 0.0:
z60 := 0.0:
V6x0 := 0.0:
V6y0 := -2.532:
V6z0 := 0.0:

x70 := 28.62:
y70 := 0.0:
z70 := 0.0:
V7x0 := 0.0:
V7y0 := -1.869:
V7z0 := 0.0:

x80 := 57.6:
y80 := 0.0:
z80 := 0.0:
V8x0 := 0.0:
V8y0 := -1.3176:
V8z0 := 0.0:

x90 := 90.3:
y90 := 0.0:
z90 := 0.0:
V9x0 := 0.0:
V9y0 := -1.05234:
V9z0 := 0.0:

xA0 := 200.0:
yA0 := 0.0:
zA0 := 0.0:
VAx0 := -0.707106:
VAy0 := 0.0:
VAz0 := 0.0:

ConditionsInitiales :=
x1(0) = x10, y1(0) = y10, z1(0) = z10, D(x1)(0) = V1x0, D(y1)(0) = V1y0, D(z1)(0) = V1z0,
x2(0) = x20, y2(0) = y20, z2(0) = z20, D(x2)(0) = V2x0, D(y2)(0) = V2y0, D(z2)(0) = V2z0,
x3(0) = x30, y3(0) = y30, z3(0) = z30, D(x3)(0) = V3x0, D(y3)(0) = V3y0, D(z3)(0) = V3z0,
x4(0) = x40, y4(0) = y40, z4(0) = z40, D(x4)(0) = V4x0, D(y4)(0) = V4y0, D(z4)(0) = V4z0,
x5(0) = x50, y5(0) = y50, z5(0) = z50, D(x5)(0) = V5x0, D(y5)(0) = V5y0, D(z5)(0) = V5z0,
x6(0) = x60, y6(0) = y60, z6(0) = z60, D(x6)(0) = V6x0, D(y6)(0) = V6y0, D(z6)(0) = V6z0,
x7(0) = x70, y7(0) = y70, z7(0) = z70, D(x7)(0) = V7x0, D(y7)(0) = V7y0, D(z7)(0) = V7z0,
x8(0) = x80, y8(0) = y80, z8(0) = z80, D(x8)(0) = V8x0, D(y8)(0) = V8y0, D(z8)(0) = V8z0,
x9(0) = x90, y9(0) = y90, z9(0) = z90, D(x9)(0) = V9x0, D(y9)(0) = V9y0, D(z9)(0) = V9z0,
xA(0) = xA0, yA(0) = yA0, zA(0) = zA0, D(xA)(0) = VAx0, D(yA)(0) = VAy0, D(zA)(0) = VAz0:

Trajectoire := dsolve({eq1x, eq1y, eq1z, eq2x, eq2y, eq2z, eq3x, eq3y, eq3z, eq4x, eq4y, eq4z, eq5x, eq5y, eq5z,eq6x, eq6y, eq6z,eq7x, eq7y, eq7z, eq8x, eq8y, eq8z, eq9x, eq9y, eq9z,eqAx, eqAy, eqAz, ConditionsInitiales}, {x1(t), y1(t), z1(t), x2(t), y2(t), z2(t), x3(t), y3(t), z3(t), x4(t), y4(t), z4(t), x5(t), y5(t), z5(t), x6(t), y6(t), z6(t), x7(t), y7(t), z7(t), x8(t), y8(t), z8(t), x9(t), y9(t), z9(t), xA(t), yA(t), zA(t)}, numeric, range = TempsInitial..TempsFinal, maxfun=0):
plots[odeplot](Trajectoire, [[x1(t), y1(t), t], [x2(t), y2(t), t], [x3(t), y3(t), t], [x4(t), y4(t), t], [x5(t), y5(t), t], [x6(t), y6(t), t], [x7(t), y7(t),t], [x8(t), y8(t), t], [x9(t), y9(t), t], [xA(t), yA(t), t]], TempsInitial..TempsFinal, numpoints = 10000,axes=boxed, scaling = constrained);

@mehdi jafari Its really long about 30 min.

@mehdi jafari done!

@mehdi jafari It's a really long code

@Carl Love I tried but I had this message:

Error, (in plots/animate) T = (24/25)*TempsInitial+(1/25)*TempsFinal .. TempsFinal is an invalid animation range

Thanks!

@Preben Alsholm Trajectoire is :

Trajectoire := dsolve({eq1x, eq1y, eq1z, eq2x, eq2y, eq2z, eq3x, eq3y, eq3z, eq4x, eq4y, eq4z, eq5x, eq5y, eq5z,eq6x, eq6y, eq6z,eq7x, eq7y, eq7z, eq8x, eq8y, eq8z, eq9x, eq9y, eq9z,eqAx, eqAy, eqAz, ConditionsInitiales}, {x1(t), y1(t), z1(t), x2(t), y2(t), z2(t), x3(t), y3(t), z3(t), x4(t), y4(t), z4(t), x5(t), y5(t), z5(t), x6(t), y6(t), z6(t), x7(t), y7(t), z7(t), x8(t), y8(t), z8(t), x9(t), y9(t), z9(t), xA(t), yA(t), zA(t)}, numeric, range = TempsInitial..TempsFinal, maxfun=0):