Problème N corps RK4 Fortran

Bonjour à tous, j'ai écris un programme fortran qui me permet d'obtenir la trajectoire de N corps en interaction gravitationnelle, j'ai choisi le soleil comme centre du repère. mais en affichant les coordonnées de la terre je me suis aperçu que c'était faux, et je ne trouve pas l'erreur, quelqu'un pourrait-il m'aider?? je met mon code plus bas ainsi que les valeurs que j'ai utilisé.

[QUOTE]
PROGRAM Gravitation

IMPLICIT NONE
REAL (KIND=8), DIMENSION(:), ALLOCATABLE :: Fx, Fy, X, Y, Vx, Vy, m
REAL (KIND=8) :: t, tfin, h
REAL (KIND=8), PARAMETER :: G=6.67259d-11, A=1.49597870691d11
INTEGER (KIND=8) :: i, N

READ *, N, t, tfin, h

ALLOCATE(X(N), Y(N), Fx(N), Fy(N), Vx(N), Vy(N), m(N))


DO i=1,N
READ *, X(i), Y(i), Vx(i), Vy(i)
END DO

DO i=1,N
READ *, m(i)
END DO

DO WHILE (t<=tfin)

CALL RK4(X, Y, Vx, Vy, m, N, t, h)

t=t+h

END DO

END PROGRAM Gravitation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE Force(Fx, Fy, X, Y, m, N, t)
IMPLICIT NONE
REAL (KIND=8), DIMENSION(1:N) :: Fx, Fy, X, Y, m, Vx, Vy
REAL (KIND=8), PARAMETER :: G=6.67259d-11
INTEGER (KIND=8) :: i, j, N
REAL (KIND=8) :: t


DO i=1,N

DO j=1,N
IF (j /= i) THEN
Fx(i) = Fx(i) - (1/G)*m(j)*m(i)*(X(i)-X(j))/(((X(j)-X(i))**2+(Y(j)-Y(i))**2)**(1.5))
Fy(i) = Fy(i) - (1/G)*m(j)*m(i)*(Y(i)-Y(j))/(((X(j)-X(i))**2+(Y(j)-Y(i))**2)**(1.5))
END IF
END DO

END DO

END SUBROUTINE Force
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE RK4(X, Y, Vx, Vy, m, N, t, h)
IMPLICIT NONE
REAL (KIND=8), DIMENSION(1:N) :: Fx, Fy, X, Y, Vx, Vy, m
REAL (KIND=8), DIMENSION(1:N) :: kx1, kx2, kx3, kx4
REAL (KIND=8), DIMENSION(1:N) :: ky1, ky2, ky3, ky4
REAL (KIND=8), PARAMETER :: G=6.67259d-11
REAL (KIND=8) :: h, t
INTEGER (KIND=8) :: p, N

kx1=0; kx2=0; kx3=0; kx4=0; ky1=0; ky2=0; ky3=0; ky4=0

Fx=0; Fy=0

CALL Force(Fx, Fy, X, Y, m, N, t)

kx1=h*Fx/m
ky1=h*Fy/m


Fx=0; Fy=0

CALL Force(Fx, Fy, X+h*0.5*Vx+0.25*h*kx1, Y+h*0.5*Vy+0.25*h*ky1, m, N, t)

kx2=h*Fx/m
ky2=h*Fy/m


Fx=0; Fy=0

CALL Force(Fx, Fy, X+h*0.5*Vx+0.25*h*kx2, Y+h*0.5*Vy+0.25*h*ky2, m, N, t)

kx3=h*Fx/m
ky3=h*Fy/m



Fx=0; Fy=0

CALL Force(Fx, Fy, X+h*Vx+0.5*h*kx3, Y+h*Vy+0.5*h*ky3, m, N, t)

kx4=h*Fx/m
ky4=h*Fy/m


X=X+h*Vx+(h/6.0d0)*(kx1+kx2+kx3)

Y=Y+h*Vy+(h/6.0d0)*(ky1+ky2+ky3)

Vx=Vx+(kx1+kx2+kx2+kx3+kx3+kx4)/6.0d0

Vy=Vy+(ky1+ky2+ky2+ky3+ky3+ky4)/6.0d0


X(N)=0; Y(N)=0; Vx(N)=0; Vy(N)=0

PRINT *, t, X(2), Y(2)

END SUBROUTINE RK4
[QUOTE]

[QUOTE]
4 0 500 0.01
-1.01723985460 -0.206180767821 0.005269822149 -0.014839159853
-0.98323814725 -0.176124724912 0.002751137188 -0.016991888551
-0.98545358338 -0.176940858196 0.002966627399 -0.017595568343
0 0 0 0
1.8
3.98600440d14
4.902798d12
1.32712440018d20
[QUOTE]