# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

def orbit(u):
    x,y,v_x,v_y = u
    r=np.hypot(x,y)
    #r= 1.521e+06
    #M,G=1.989e+30,6.7e-11
    M,G=20,110
    f=G*M/r**3
    return np.array([v_x,v_y,-f*x,-f*y])


def RK4(f,u,dt):
    k1=f(u)*dt
    k2=f(u+0.5*k1)*dt
    k3=f(u+0.5*k2)*dt
    k4=f(u+k3)*dt
    return u+(k1+2*k2+2*k3+k4)/6

    
def RK4_int(f,y0,tspan):
    y=np.zeros([len(tspan),len(y0)])
    y[0,:] =y0
    for k in range (1,len(tspan)):
        y[k,:] = RK4(f,y[k-1],tspan[k]-tspan[k-1])
    return y

dt=0.1
t = np.arange(0,10,dt)
y0=np.array([10, 0.0, 10, 10])

sol_rk4=RK4_int(orbit,y0,t)
x,y,v_x,v_y = sol_rk4.T
plt.grid()
plt.plot(x,y)
plt.show()

[((31, 4, 31, 22), 'numpy.arange', 'np.arange', ({(31, 14, 31, 15): '0', (31, 16, 31, 18): '10', (31, 19, 31, 21): 'dt'}, {}), '(0, 10, dt)', True, 'import numpy as np\n'), ((32, 3, 32, 30), 'numpy.array', 'np.array', ({(32, 12, 32, 29): '[10, 0.0, 10, 10]'}, {}), '([10, 0.0, 10, 10])', True, 'import numpy as np\n'), ((36, 0, 36, 10), 'matplotlib.pyplot.grid', 'plt.grid', ({}, {}), '()', True, 'import matplotlib.pyplot as plt\n'), ((37, 0, 37, 13), 'matplotlib.pyplot.plot', 'plt.plot', ({(37, 9, 37, 10): 'x', (37, 11, 37, 12): 'y'}, {}), '(x, y)', True, 'import matplotlib.pyplot as plt\n'), ((38, 0, 38, 10), 'matplotlib.pyplot.show', 'plt.show', ({}, {}), '()', True, 'import matplotlib.pyplot as plt\n'), ((7, 6, 7, 19), 'numpy.hypot', 'np.hypot', ({(7, 15, 7, 16): 'x', (7, 17, 7, 18): 'y'}, {}), '(x, y)', True, 'import numpy as np\n'), ((12, 11, 12, 40), 'numpy.array', 'np.array', ({(12, 20, 12, 39): '[v_x, v_y, -f * x, -f * y]'}, {}), '([v_x, v_y, -f * x, -f * y])', True, 'import numpy as np\n')]