RK4 Second order ODE debugging

So I have been tasked with solving a set of ODE's using the RK4 method and while my code runs I don't believe I am getting the right values. My hunch is that my error resides in the way I am initializing my k values and the way I create my intermediate equations.

here is my code, (I know its long, not too great with condensing code in matlab just yet, sorry)

ive been told the final x value should be about -.92359 and the final y value should be .60566

ill also just attach it as an .mfile if true % code clear; clc; x(1) = 1.2; y(1) = 0; vx(1) = 0; vy(1) = -1.0493571; dx(1) = vx(1); dy(1) = vy(1); t0 = 0; tf = 10; f = 0; N = 1000; dt = (tf-t0)/N; eps = 0.01; mu1 = 1/82.45; mu2 = 1-mu1; T = t0:dt:tf;

r1 = (((x((((x(i)+mu1)^2 + y(i)^2)^.5);i)+mu1)^2 + y(i)^2)^.5); r2 = (((x(i)-mu2)^2 +y(i)^2)^.5);

%F(x,y,dx,dy) = @(x,y,dx,dy) [2*dy - f*dx + x - (mu2*(x+mu1))/(((x+mu1)^2 +y^2)^1.5) - (mu1*(x-mu2))/(((x-mu2)^2 +y^2)^1.5), -2*dx - f*dy + y - (-mu2*y)/(((x+mu1)^2 +y^2)^1.5) - (-mu1*y)/(((x-mu2)^2 +y^2)^1.5)];

%X(x,y,dx,dy) = @(x,y,dx,dy) 2*dy - f*dx + x - (mu2*(x+mu1))/(((x+mu1)^2 +y^2)^1.5) - (mu1*(x-mu2))/(((x-mu2)^2 +y^2)^1.5);

%Y(x,y,dx,dy) = @(x,y,dx,dy) -2*dx - f*dy + y - (-mu2*y)/(((x+mu1)^2 +y^2)^1.5) - (-mu1*y)/(((x-mu2)^2 +y^2)^1.5);

for  i=1:(length(T)-1);

 %while (r1(i+1)-r1(i))/(r1(i+1)) < eps   
  %initial time step t=0 set the initial values for Fs  
  F1(1) = x(1);
  F2(1) = dx(1);
  F3(1) = y(1);
  F4(1) = dy(1);

x(i)  =  F1(i);
dx(i) =  F2(i);
y(i)  =  F3(i);
dy(i) =  F4(i);

%begin creating K1s to create intermediate F values
k11 = F2(i);
k12 = 2*F4(i) - f*F2(i) +F1(i) - mu2*(F1(i)+mu1)/(((F1(i)+mu1)^2 + F2(i)^2)^1.5) - mu1*(F1(i)-mu2)/(((F1(i)-mu2)^2 + F3(i)^2)^1.5);
k13 = F4(i);
k14 = -2*F2(i) - f*F4(i) + F3(i) -mu2*F3(i)/(((F1(i)+mu1)^2 + F3(i)^2)^1.5) - mu1*F3(i)/(((F1(i)-mu2)^2 + F3(i)^2)^1.5);

%initialize values for intermediate F equations (intermediate slopes)
F1tmp(i) = F1(i) + (dt/2)*k11;
F2tmp(i) = F2(i) + (dt/2)*k12; 
F3tmp(i) = F3(i) + (dt/2)*k13; 
F4tmp(i) = F4(i) + (dt/2)*k14;

%create K2s for next intermediate F values
k21 = F2tmp(i);
k22 = 2*F4tmp(i) - f*F2tmp(i) +F1tmp(i) - mu2*(F1tmp(i)+mu1)/(((F1tmp(i)+mu1)^2 + F2tmp(i)^2)^1.5) - mu1*(F1tmp(i)-mu2)/(((F1tmp(i)-mu2)^2 + F3tmp(i)^2)^1.5);
k23 = F4tmp(i);
k24 = -2*F2tmp(i) - f*F4tmp(i) + F3tmp(i) -mu2*F3tmp(i)/(((F1tmp(i)+mu1)^2 + F3tmp(i)^2)^1.5) - mu1*F3tmp(i)/(((F1tmp(i)-mu2)^2 + F3tmp(i)^2)^1.5);

%refresh intermediate F values for next set of Ks
F1tmp(i) = F1tmp(i) + (dt/2)*k21;
F2tmp(i) = F2tmp(i) + (dt/2)*k22; 
F3tmp(i) = F3tmp(i) + (dt/2)*k23; 
F4tmp(i) = F4tmp(i) + (dt/2)*k24;

%create K3s for next set of intermediate F values
k31 = F2tmp(i);
k32 = 2*F4tmp(i) - f*F2tmp(i) +F1tmp(i) - mu2*(F1tmp(i)+mu1)/(((F1tmp(i)+mu1)^2 + F2tmp(i)^2)^1.5) - mu1*(F1tmp(i)-mu2)/(((F1tmp(i)-mu2)^2 + F3tmp(i)^2)^1.5);
k33 = F4tmp(i);
k34 = -2*F2tmp(i) - f*F4tmp(i) + F3tmp(i) -mu2*F3tmp(i)/(((F1tmp(i)+mu1)^2 + F3tmp(i)^2)^1.5) - mu1*F3tmp(i)/(((F1tmp(i)-mu2)^2 + F3tmp(i)^2)^1.5);

%refresh intermediate F values for final K values
F1tmp(i) = F1tmp(i) + (dt/2)*k31;
F2tmp(i) = F2tmp(i) + (dt/2)*k32; 
F3tmp(i) = F3tmp(i) + (dt/2)*k33; 
F4tmp(i) = F4tmp(i) + (dt/2)*k34;

%create final K values to use in calculating F(i+1)
k41 = F2tmp(i);
k42 = 2*F4tmp(i) - f*F2tmp(i) +F1tmp(i) - mu2*(F1tmp(i)+mu1)/(((F1tmp(i)+mu1)^2 + F2tmp(i)^2)^1.5) - mu1*(F1tmp(i)-mu2)/(((F1tmp(i)-mu2)^2 + F3tmp(i)^2)^1.5);
k43 = F4tmp(i);
k44 = -2*F2tmp(i) - f*F4tmp(i) + F3tmp(i) -mu2*F3tmp(i)/(((F1tmp(i)+mu1)^2 + F3tmp(i)^2)^1.5) - mu1*F3tmp(i)/(((F1tmp(i)-mu2)^2 + F3tmp(i)^2)^1.5);

%calculate F values at next time step
F1(i+1) = F1(i) + (dt/6)*(k11 + 2*k21 + 2*k31 + k41);
F2(i+1) = F2(i) + (dt/6)*(k12 + 2*k22 + 2*k32 + k42);  
F3(i+1) = F3(i) + (dt/6)*(k13 + 2*k23 + 2*k33 + k43);  
F4(i+1) = F4(i) + (dt/6)*(k14 + 2*k24 + 2*k34 + k44);

end end