clear; clc; close all;
Alpha=1; dt=0.025;
x1=2.5; y1=0;
vx1=0; vy1=0.25;
Emax=-1e18;
Emin=1e18;
r=sqrt(x1^2+y1^2); 
E1=(vx1^2+vy1^2)/4 ...
- Alpha/r; 
hl=plot(x1,y1,'o');
a=Alpha/(2*abs(E1));
b=r*vy1/sqrt(2*abs(E1)); 
axis([-2*a 2*a -1.2*b 1.2*b]); 
set(hl,'Color','g');
ha=gca;
set(ha,'XTick',[-a 0 a],'YTick',[-b 0 b]);
grid on; 
pause 
while 1
x2=x1+vx1*dt; 
y2=y1+vy1*dt; 
r=sqrt(x2^2+y2^2); 
A = Alpha/r^2; 
ax=-A*x2/r; 
ay=-A*y2/r; 
vx2=vx1+ax*dt; 
vy2=vy1+ay*dt; 
E2=((vx2+vx1)^2+(vy2+vy1)^2)/8 ...
-Alpha/r; 
if E2> Emax 
Emax=E2; 
end;
if E2< Emin
Emin=E2;
end;
set(hl,'XData',[x1 x2],'YData',[y1 y2]);
pause(0.01)
x1=x2;
y1=y2;
vx1=vx2;
vy1=vy2;
end;