Моделирование орбитального движения вокруг тела в Matlab

Я работаю над очень простой программой для имитации орбитального движения тела вокруг другого тела, как спутник вокруг Земли. Я следовал рекомендациям и уравнениям, доступным в книгах и Интернете, но объект, похоже, вообще не вращается по орбите. Пожалуйста, помогите мне, если кто-то может. Заранее спасибо. Код ниже

v0=80;
theta=45*pi/180;

vx(1)=v0*cos(theta);
vy(1)=v0*sin(theta);
px(1)=0;
py(1)=0;

mass=100; %kg
cmass=400; % mass of the body at 400,500
    ax=0;
%
g=9.8 ;%m/s^2
ay=-g;

p2x=400; % x co-ordinate of the stationary body
p2y=500; % y co-ordinate of the stationary body

G=6.674*10^-11; % the Gravitational Constant
figure(1)
plotsize=800;
i=1;
dt=.1;
t=0;
while(t<20)

%a2x=a2x-0.10;
%a2y=a2y+0.50;
r=sqrt((p2x-px(i))^2+(p2y-py(i))^2); % distance between the two bodies
F=((G*mass*cmass)/r^2);   % force by formula f=(G*m1*m2)/r^2
a=-(1/cmass^2)*F;    % acceleraion  a=1/m1^2*F

vx(i+1)=vx(i)+(ax)*dt;
vy(i+1)=vy(i)+(ay)*dt;

px(i+1)=px(i)+vx(i)*dt;
py(i+1)=py(i)+vy(i)*dt;

hold off
plot(px(i+1),py(i+1),'o','MarkerSize',15)

hold on
plot(px,py,'r')
plot(p2x,p2y);
axis([0 plotsize 0 plotsize])
pause(.1) %pause for graphics
i=i+1;
t=t+dt;

конец