Make Fluxion with the gratiational potential energy beween body i and j

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