If you drop `N`, you will see that the result is a very long and complex expression, full of high-degree `Root` objects. I would say that the result isn't really a bug, but shows an accumulation of numerical errors, which you can avoid by using arbitrary-precision calculations (set by the second argument of `N`), for example,

````
Maximize[{f[x1, x2, x3], x1^2 + x2^2 + x3^2 <= 1}, {x1, x2, x3}] // N[#, 5]&
(* {0.057092, {x1 -> -0.87364, x2 -> -0.026513, x3 -> -0.48585}} *)
````

Alternatively, you can use numerical maximization with [`NMaximize`](https://reference.wolfram.com/language/ref/NMaximize.html):

````
NMaximize[{f[x1, x2, x3], x1^2 + x2^2 + x3^2 <= 1}, {x1, x2, x3}]
(* {0.0570927, {x1 -> -0.873678, x2 -> -0.0266629, x3 -> -0.485778}} *)