x1 = -1;
y1 = 0;
vx1 = 0;
vy1 = -0.2;

x2 = 1;
y2 = 0;
vx2 = 0;
vy2 = 0.1;


function dx = f(X, t)
m1 = 1;
m2 = 2;

x1 = X(1);
y1 = X(2);
vx1 = X(3);
vy1 = X(4);
x2 = X(5);
y2 = X(6);
vx2 = X(7);
vy2 = X(8);

Dx = x2 - x1;
Dy = y2 - y1;
r  = sqrt(Dx^2 + Dy^2);
wsp = m1*m2/(r^3);
Fx = wsp * Dx;
Fy = wsp * Dy;

dx(1) = vx1;
dx(2) = vy1;
dx(3) = Fx/m1;
dx(4) = Fy/m1;

dx(5) = vx2;
dx(6) = vy2;
dx(7) = -Fx/m2;
dx(8) = -Fy/m2;

endfunction

X0 = [x1, y1, vx1, vy1, x2, y2, vx2, vy2];
t = 0:0.1:6;

y = lsode("f", X0, t);
plot(y(1,1), y(1,2), "rx");
plot(y(1,5), y(1,6), "bo");
axis ([-1,1, -0.4, 0.4]);
hold on ;
for i = 2:length(t)
  plot(y(i,1), y(i,2), "rx");
  plot(y(i,5), y(i,6), "bo");
  pause (0.05);
endfor

hold off;



