lagrange_points should support vectorized inputs

For example, to compute for arrays of position and velocity (and time).

Failing example:
```python
gd.cluster.lagrange_points(
    pot, orbit_q, orbit_p, u.Quantity(1e5, "Msun"), t_grid
)
```
where `orbit_q.shape == orbit_p.shape == (N, 3)` and `t_grid.shape == (N,)`