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

x2(t):=cos(t):
y2(t):=sin(t):

Position := [x(t), y(t), z(t)]:
#Soleil:=[1,2,0]:
Soleil:=[x2(t),y2(t),0]:

r := sqrt(x(t)^2 + y(t)^2 + z(t)^2 + 0.001):
r2:= sqrt((x2(t)-x(t))^2 + (y2(t)-y(t))^2 + (0)^2 + 0.001):

Vitesse := [(D(x))(t), (D(y))(t), (D(z))(t)]:

CouplageGravitation := 1.0:
CouplageGravitation2 := 0.005:

ForceGravitationnelle := -CouplageGravitation*Position/r^3:
ForceGravitationnelle2:= -CouplageGravitation2*Soleil/(r2)^3:


Force := ForceGravitationnelle + ForceGravitationnelle2:

Fx := innerprod([1, 0, 0], Force);
Fy := innerprod([0, 1, 0], Force);
Fz := innerprod([0, 0, 1], Force);

TempsInitial := 0:
TempsFinal := 250:

eq1 := (D(D(x)))(t) = Fx;
eq2 := (D(D(y)))(t) = Fy;
eq3 := (D(D(z)))(t) = Fz;
  /
- |         /      2                       2         2                       2
  \1.0 x(t) \cos(t)  - 2 cos(t) x(t) + x(t)  + sin(t)  - 2 sin(t) y(t) + y(t)

                                                                (3/2)\//
          \                      /    2       2       2        \     | |
   + 0.001/^(3/2) + 0.005 cos(t) \x(t)  + y(t)  + z(t)  + 0.001/     / \

                                 (3/2)
  /    2       2       2        \      /      2                       2
  \x(t)  + y(t)  + z(t)  + 0.001/      \cos(t)  - 2 cos(t) x(t) + x(t)

                                                   \
           2                       2        \      |
   + sin(t)  - 2 sin(t) y(t) + y(t)  + 0.001/^(3/2)/
  /
- |         /      2                       2         2                       2
  \1.0 y(t) \cos(t)  - 2 cos(t) x(t) + x(t)  + sin(t)  - 2 sin(t) y(t) + y(t)

                                                                (3/2)\//
          \                      /    2       2       2        \     | |
   + 0.001/^(3/2) + 0.005 sin(t) \x(t)  + y(t)  + z(t)  + 0.001/     / \

                                 (3/2)
  /    2       2       2        \      /      2                       2
  \x(t)  + y(t)  + z(t)  + 0.001/      \cos(t)  - 2 cos(t) x(t) + x(t)

                                                   \
           2                       2        \      |
   + sin(t)  - 2 sin(t) y(t) + y(t)  + 0.001/^(3/2)/
                                   1.0 z(t)
                   - ------------------------------------
                                                    (3/2)
                     /    2       2       2        \
                     \x(t)  + y(t)  + z(t)  + 0.001/
                      /
   @@(D, 2)(x)(t) = - |         /      2                       2         2
                      \1.0 x(t) \cos(t)  - 2 cos(t) x(t) + x(t)  + sin(t)

                            2        \
      - 2 sin(t) y(t) + y(t)  + 0.001/^(3/2)

                                                    (3/2)\//
                     /    2       2       2        \     | |
      + 0.005 cos(t) \x(t)  + y(t)  + z(t)  + 0.001/     / \

                                    (3/2)
     /    2       2       2        \      /      2                       2
     \x(t)  + y(t)  + z(t)  + 0.001/      \cos(t)  - 2 cos(t) x(t) + x(t)

                                                      \
              2                       2        \      |
      + sin(t)  - 2 sin(t) y(t) + y(t)  + 0.001/^(3/2)/
                      /
   @@(D, 2)(y)(t) = - |         /      2                       2         2
                      \1.0 y(t) \cos(t)  - 2 cos(t) x(t) + x(t)  + sin(t)

                            2        \
      - 2 sin(t) y(t) + y(t)  + 0.001/^(3/2)

                                                    (3/2)\//
                     /    2       2       2        \     | |
      + 0.005 sin(t) \x(t)  + y(t)  + z(t)  + 0.001/     / \

                                    (3/2)
     /    2       2       2        \      /      2                       2
     \x(t)  + y(t)  + z(t)  + 0.001/      \cos(t)  - 2 cos(t) x(t) + x(t)

                                                      \
              2                       2        \      |
      + sin(t)  - 2 sin(t) y(t) + y(t)  + 0.001/^(3/2)/
                                            1.0 z(t)
           @@(D, 2)(z)(t) = - ------------------------------------
                                                             (3/2)
                              /    2       2       2        \
                              \x(t)  + y(t)  + z(t)  + 0.001/
x0 := 2.0:
y0 := 0.1:
z0 := 1.0:
Vx0 := -0.5:
Vy0 := 0.0:
Vz0 := 0.0:

ConditionsInitiales := x(0) = x0, y(0) = y0, z(0) = z0, D(x)(0) = Vx0, D(y)(0) = Vy0, D(z)(0) = Vz0 ;
Trajectoire := dsolve({eq1, eq2, eq3, ConditionsInitiales}, {x(t), y(t), z(t)}, numeric, range = TempsInitial..TempsFinal, maxfun=0);
x(0) = 2.0, y(0) = 0.1, z(0) = 1.0, D(x)(0) = -0.5, D(y)(0) = 0., D(z)(0) = 0.
proc(x_rkf45)  ...  end;
plots[odeplot](Trajectoire, [x(t), y(t), z(t)], TempsInitial..TempsFinal, numpoints = 1000, axes = boxed, scaling = constrained);