Solve with Python please!

Consider two masses m_1, m_2 which are at time t in the positions x_1 (t), x_2 (t). According to Newton's law of gravitation the mass m_2 attracts the mass m_1 with the force F_1 (x_1 (t), x_2 (t) = gamma m_1 m_2/||x_2 (t) - x_1 (t)||_2^3 (x_2 (t) - x_1 (t)), where gamma is the gravity constant. The attraction force of m_1 on m_2 is F_2 = -F_1. Let v_k (t) = x_k (t) the velocity of the mass k at time t. Then m_k v_k (t) = F_k (t), k = 1,2 (force - mass times acceleration). The motion of the masses is described by the system of ordinary differential equations d/dt [x_1 (t) x_2 (t) v_1 (t) v_2 (t)] = [v_1 (t) v_2 (t) f_1 (x_1 (t), x_2 (t))/m_1 F_2 (x_1 (t), x_2 (t))/m_2] Write a Matlab or Python function twomasses(m1,m2,x1,x2,p,h) which solves the system (*) using a classical Runge-Kutta method and generate an animation showing the motion of the masses. Assume all vectors are 2-dimensional (plane motion). The parameter h is the time step size, and the parameter p determines the initial velocities: v_1 (t_0) = [0 p/m_1], v_2 (t_0) = [0 -p/m_2] Assume the gravitational constant is gamma = 1