input: Write a function to simulate the motion of two masses under gravitational attraction and return their positions at each time step.
  assert np.allclose(simulate_two_body_motion(1, 1, 0, 0, 1, 0, 3, 0, 0, 0.5, 1000, 1e-3)[0][-1], [-0.49949056,  0.        ], rtol=1e-5)
  assert np.allclose(simulate_two_body_motion(2, 2, 0, 0, 1, 0, 3, 0, 0, 0.5, 1000, 1e-3)[1][-1], [ 3.49949056,  0.        ], rtol=1e-5)
  assert np.allclose(simulate_two_body_motion(1, 1, 0, 0, 1, 0, 3, 1, 0, 0.5, 1000, 1e-3)[0][100], [ 0.099995,  0.        ], rtol=1e-5)
output: import numpy as np
  
  def simulate_two_body_motion(m1, m2, x1, y1, vx1, vy1, x2, y2, vx2, vy2, G, steps, dt):
      positions1 = np.zeros((steps, 2))
      positions2 = np.zeros((steps, 2))
      positions1[0] = [x1, y1]
      positions2[0] = [x2, y2]
      
      for i in range(1, steps):
          r12 = np.array([positions2[i-1, 0] - positions1[i-1, 0], positions2[i-1, 1] - positions1[i-1, 1]])
          r = np.sqrt(r12[0]**2 + r12[1]**2)
          F = G * m1 * m2 * r12 / r**3
          
          ax1 = F[0] / m1
          ay1 = F[1] / m1
          ax2 = -F[0] / m2
          ay2 = -F[1] / m2
          
          vx1 += ax1 * dt
          vy1 += ay1 * dt
          vx2 += ax2 * dt
          vy2 += ay2 * dt
          
          x1 += vx1 * dt
          y1 += vy1 * dt
          x2 += vx2 * dt
          y2 += vy2 * dt
          
          positions1[i] = [x1, y1]
          positions2[i] = [x2, y2]
      
      return positions1, positions2