I have a multivariate (in practice, 2β5 dimensions) function approximated with Chebyshev polynomials. I would like to find *all* the local extrema of the approximation. This is equivalent to solving the system obtained from partial derivatives, but there may be a more direct method (also, dealing with extrema at endpoints, these are bounded on [-1,1]^n.

Specifically, I have the approximation available as either weighted sums of a Chebyshev basis, or if preferable vectors of

```
(i j k ...) => c
```

where c would be the coefficient assigned to x_1^i x_2^j x_3^k \dots.

Is this a tractable problem? I only know a method for the unviariate case (Boyd 2013).

If yes, what would be the recommended way in Julia?