The band extrema may be at an arbitrary point in the Brillouin zone, so determining their position can become computationally demanding. The following strategy is something that is sometimes used to locate Weyl points (band crossings) in a band structure, which I think should also be useful for locating band extrema. The strategy would be:
- Sample the full Brillouin zone with a uniform $\mathbf{k}$-point grid.
- Locate the band extrema in your $\mathbf{k}$-point grid.
- Sample a patch of the Brillouin zone around the band extrema of your previous grid.
- Iterate until convergence of the location of the band extrema.
There are still some choices to be made (density of sampling grids, size of patches), and I imagine this will require some trial-and-error as ideal values will depend on the specific band structure (effective masses, bandwidth).
The approach does sound computationally demanding, but you should be able to diagonalize the Hamiltonian over the required $\mathbf{k}$-point grids in a non-self-consistent manner, which should accelerate the calculations.