I have a problem that I need to simulate where there are four particles located at the corners of a square room with known lenght a. They are starting to move simultaneously and are supposed to move towards the closest neighbour until one of them reaches some other. They all have constant velocity V. I am supposed to use scipy.integrate.RK45 function.

I am not sure how to describe mathematically the fact that they are all moving towards the closest neightbour. I suppose the function should look something like this (for one particles for now):

def motion_eq(t, y):
        ???

r = np.array([0, 0], 'float')
v = np.array([0, 1], 'float')

dt = 0.1
a = 2

motion_solver = scipy.integrate.RK45(motion_eq, 0, np.hstack([r, v]),
    t_bound = np.inf, first_step = dt, max_step = dt)

particle, *_ = plt.plot(*r.T, 'o')
plt.gca().set_aspect(1)
plt.xlim([-a, a])
plt.ylim([-a, a])
def update():
    motion_solver.step()
    r = motion_solver.y[0:2]
    particle.set_data(*r.T)
    plt.draw()
timer = plt.gcf().canvas.new_timer(interval = 50)
timer.add_callback(update)
timer.start()

plt.show()

Could somebody please give a hint of how to do this? Thanks in advance.