if true
f= @(R_moon) (-G*(moon_mass+earth_mass)*R_moon)/norm(R_moon)^3;
R1=moon_position_vector(2,:);
V1=data_moon(2,8:10);
t=1; tfin=24*3600;
h=1; hmin=h/64; hmax=64*h;
j=1;
epsilon=0.0000001;
max1=200;
while (t<=tfin)
k1 = h*f(R1);
k2 = h*f(R1+(1/4)*k1);
k3 = h*f(R1+(3/32)*k1+(9/32)*k2);
k4 = h*f(R1+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3);
k5 = h*f(R1+(438/216)*k1-8*k2+(3680/513)*k3-(845/4104)*k4);
k6 = h*f(R1-(8/27)*k1+2*k2-(3544/2565)*k3+(1859/4104)*k4-(11/40)*k5);
err=max(abs(((1/360)*k1-(128/4275)*k3-(2197/75240)*k4+(1/50)*k5+(2/55)*k6)));
Vnew= V1+(25/216)*k1+(1408/2565)*k3+(2197/4101)*k4-(1/5)*k5;
Vnew2 = V1 + 16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - 9*k5/50 + 2*k6/55;
R=max(abs(Vnew-Vnew2));
delta=0.84*(epsilon*h/R)^(0.25);
if ((err<epsilon)|(h<2*hmin))
V(j,:)=Vnew2;
R1=R1+Vnew*h;
Rnew(j,:)=R1;
end
t=t+h;
time(j) = t;
j=j+1;
if ((delta<0.75)&(h>2*hmin)), h = h/2; end
if ((delta>1.50)&(2*h<hmax)), h = 2*h; end
end
end
