Orbital.sagews

G = 5 # These are made up numbers
M = 5
delT = 0.01

x0 = 5
y0 = 0
r0 = sqrt(x0^2+y0^2)

vx0 = 0
vy0 = 1

ax0 = G*M*x0/r0^3
ay0 = G*M*y0/r0^3

x = [x0]
y = [y0]
vx = [vx0]
vy = [vy0]
ax = [ax0]
ay = [ay0]

for i in range(5000):
    axOld = ax[-1]
    ayOld = ay[-1]
    vxOld = vx[-1]
    vyOld = vy[-1]
    
    vxNew = vxOld - delT*axOld
    vyNew = vyOld - delT*ayOld

    xOld = x[-1]
    yOld = y[-1]

    xNew = xOld + delT*vxNew
    yNew = yOld + delT*vyNew
    
    rNew = sqrt(xNew^2+yNew^2)
    axNew = G*M*xNew/rNew^3
    ayNew = G*M*yNew/rNew^3
    
    x.append(xNew)
    y.append(yNew)
    vx.append(vxNew)
    vy.append(vyNew)
    ax.append(axNew)
    ay.append(ayNew)
    
list_plot(zip(x,y))