docs: def U_ij(q_vars: List[fl.Var], mass: np.ndarray, i: int, j: int):
      """Make Fluxion with the gratiational potential energy beween body i and j"""
      # Check that the lengths are consistent
      assert len(q_vars) == 3 * len(mass)
      
      # Masses of the bodies i and j
      mi = mass[i]
      mj = mass[j]
      
      # Gravitational potential is -G * m1 * m2 / r
      U = -(G * mi * mj) / flux_r(q_vars, i, j)
      return U