// (x1,y1) = (-0.97000436, 0.24308753), (x2,y2) = (-x1, -y1), (x3,y3) = (0,0) 
// (vx1,vy1) = (vx2, vy2) = -(vx3, vy3)/2; where (vx3,vy3) = (0.93240737, 0.86473146) 
// m=1; G=1

function dz=F(t,z)
    dx12=z(1)-z(5);dx13=z(1)-z(9);dx23=z(5)-z(9)
    dy12=z(3)-z(7);dy13=z(3)-z(11);dy23=z(7)-z(11)
    d12=(dx12^2+dy12^2)^1.5;d13=(dx13^2+dy13^2)^1.5
    d23=(dx23^2+dy23^2)^1.5
    for i=0:5
        dz(2*i+1)=z(2*i+2)
    end
    dz(2)=-dx12/d12-dx13/d13
    dz(4)=-dy12/d12-dy13/d13
    dz(6)=dx12/d12-dx23/d23
    dz(8)=dy12/d12-dy23/d23
    dz(10)=dx13/d13+dx23/d23
    dz(12)=dy13/d13+dy23/d23
endfunction


x1=-0.97000436;x2=-x1
y1=0.24308753;y2=-y1
x3=0;y3=0
vx3=0.93240737;vy3=0.86473146;
//vx3=0.6;vy3=0.86473146
vx2=-vx3/2;vy2 =-vy3/2
vx1=vx2;vy1=vy2


cond_in=[x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3]'

n=10000
h=(10-0)/n
sol=zeros(12,n+1)
t=0
sol(:,1)=cond_in
for i=1:n
    sol(:,i+1)=RK4(sol(:,i),t,t+h,F)
    t=t+h
end

scf(0)
plot(sol(1,:),sol(3,:),'.r')
plot(sol(5,:),sol(7,:),'.g')
plot(sol(9,:),sol(11,:),'.b')
xgrid
