clear;  clc;

%    selftest03:  2 masses under gravitation
%    our state vector will be (x1,y1,vx1,vy1)


G=6.67e-11;
M=1;  G=1;
yinitial=[ -1  0  0  1.3 ];
tstart=0; tend=100; dt=0.001;

k=0;
y=yinitial;
for time=tstart:dt:tend
    k=k+1;
    xx1(k)=y(1);
    yy1(k)=y(2);
    yprime=RHS(y,G,M);
    ytilde = y + dt*yprime;
    yprimetilde=RHS(ytilde,G,M);
    ytildetilde=y + dt*yprimetilde;
    y = 0.5*( ytilde + ytildetilde );
end

plot(xx1,yy1);
fprintf('done after %d steps\n',k);;

function dydt=RHS(y,G,M)
   R=sqrt( y(1)^2 + y(2)^2 );  R3i=1/(R^3);
   dydt(1)=y(3);
   dydt(2)=y(4);
   dydt(3)=-G*M*y(1)*R3i;
   dydt(4)=-G*M*y(2)*R3i;
end
