function planet(eps,N)
% integrate equations of planetary motion
% by a 2nd order leap-frog scheme
[a,b,c]=sphere(100);
r(1,1) = -.5;r(1,2)=0;v0=[0,-1.63];
ep2=eps/2;
v(1,:)=v0+ep2*accel(r(1,1),r(1,2));
for i=1:N
  r(i+1,:) = r(i,:) + eps*v(i,:);
  v(i+1,:) = v(i,:) + eps*accel(r(i+1,1),r(i+1,2));
end
figure; axis([-.55, 1.05 , -.8, .8]); hold on
contour(.1*a,.1*b,.1*c)
plot(r(:,1),r(:,2),'-b',...
    'LineWidth',2)
plot(r(1:14:N,1),r(1:14:N,2),'ro','MarkerEdgeColor','k',...
    'MarkerFaceColor',[.49 1 .63],...
    'MarkerSize',12)
s='planetary motion around the sun';title(s);hold off
function a = accel(x,y)
   r3 = (x^2+y^2)^(3/2);
   a(1,1) = -x/r3; 
   a(1,2) = -y/r3; 