def findenergy(x, y, vx, vy, mu):
    """Finds the energy at a point in the 4D phase-space."""
    r1 = ((x + mu)**2 + y**2) ** 0.5 
    r2 = ((x - (1 - mu))**2 + y**2) ** 0.5
    return 0.5 * (vx**2 + vy**2 - x**2 - y**2) - (1 - mu)/r1 - mu/r2 - 0.5 * (1 - mu) * mu