Runge kutta de orden cuarto "problema de 2 cuerpos"

Program Mecanica_EULER

Implicit none

Real(8) :: Dt,vx,vy,vz,x,y,z,a,b,mu,cociente,e0,s,u,alpha,p,t
REAL(8) :: ax, ay, az, kx1, kx2, kx3, kx4, qx1, qx2, qx3, qx4, ky1, ky2, ky3, ky4,qy1, qy2, qy3, qy4
REAL(8) :: kz1, kz2, kz3, kz4, qz1, qz2, qz3, qz4
Integer(8) :: i , N

!las unidades están en el sistema astronomico, no MKS! NO CONFUNDIRalpha=8.904307867d0*10.d0**(-10) !Gm1m2
mu=3.000660266d0*10.d0**(-6) ! masa reducida
s=1.884459221d0 !suma de cuadrados de condiciones iniciales de posicion a la tres medios
x=-1.799915723d0
y=0d0
z=0d0
vx=0d0vy=-5.74085227d0*10d0**(-3)
vz=0d0
e0=-0.000000000445259975343064

ax=(-alpha)*x/(mu*s)
ay=(-alpha)*y/(mu*s)
az=(-alpha)*z/(mu*s)
!e0=0.5d0*mu*(vx**2.+vy**2.+vz**2.)-alpha/(s**(1d0/3d0))
!stopwrite(*,*)"Ingrese el valor del tiempo inicial, en días:"
read(*,*)a
write(*,*)"Ingrese el valor del tiempo final, en días:"
read(*,*)b
write(*,*)"Ingrese el paso de tiempo, al que llamaremos Dt:"read(*,*)Dt

N=(b-a)/Dt

open(25,file="aceleracionprueba.dat")
open(26,file="velocidad.dat")
open(27,file="posicion.dat")
open(28,file="energia.dat")

!WRITE(*,'(3E12.5)') x, y, z!WRITE(*,'(3E12.5)') vx, vy, vz

!CALL aceleracion(x,y,z,ax,ay,az)
!WRITE(*,*) ax
!STOP

do i=0,N

!k para posición, q para velocidades

kx1=Dt*vx
ky1=Dt*vy
kz1=Dt*vz

qx1=Dt*axqy1=Dt*ay
qz1=Dt*az


CALL aceleracion(x+kx1/2.,y+ky1/2.,z+kz1/2.,ax,ay,az)

!vx=vx+ax*Dt
!vy=vy+ay*Dt
!vz=vz+az*Dt

kx2=Dt*(vx+qx1/2.)
ky2=Dt*(vy+qy1/2.)kz2=Dt*(vz+qy1/2.)

qx2=Dt*ax
qy2=Dt*ay
qz2=Dt*az

!write(*,*) vz
!stop

CALL aceleracion(x+kx2/2.,y+ky2/2.,z+kz2/2.,ax,ay,az)

!vx=vx+ax*Dt
!vy=vy+ay*Dt
!vz=vz+az*Dtkx3=Dt*(vx+qx2/2.)
ky3=Dt*(vy+qy2/2.)
kz3=Dt*(vz+qz2/2.)

qx3=Dt*ax
qy3=Dt*ay
qz3=Dt*az

CALL aceleracion(x+kx3,y+ky3,z+kz3,ax,ay,az)

!vx=vx+ax*Dt
!vy=vy+ay*Dt...