function twobody
tspan = [0;10];
y0 = [1;0;0;1];

options = odeset('RelTol',1e-6,'AbsTol',1e-6);
%options = [];
[t,y] = ode45(@f,tspan,y0,options);

plot(y(:,1),y(:,3),'.-b',0,0,'*r');

% ----------------------------------------------
function ydot = f(t,y)
%G = 6.673E-11;
%M = 2.0E30;
%g = -(y(1)^2+y(3)^2)^1.5;
g = -(y(1)^2 + y(3)^2)^-1.5;
%g = -1;
ydot = [y(2);g*y(1);y(4);g*y(3)];