First expand the function:
$$f(x, y, z) = \left(x -\frac{1}{3}\right)^2 + \left(y - \frac{1}{3}\right)^2 + \left(z - \frac{1}{3}\right)^2$$
$$f(x, y, z) = x^2 - \frac 23 x + \frac 19 + y^2 -\frac 23 y + \frac 19 + z^2 - \frac 23 z + \frac 19$$
$$f(x, y, z) = x^2 + y^2 + z^2 - \frac 23 x - \frac 23 y - \frac 23z + \frac 13$$
To minimize/maximize a function with a given constain it's best to apply Lagrange multiplier. So after applying it the function will look like this:
$$F(x, y, z, \lambda) = x^2 + y^2 + z^2 - \frac 23 x - \frac 23 y - \frac 23z + \frac 13 + \lambda (x+y+z-1)$$
Now we take partial derivatives and we set them to $0$.
$$F_x = 2x - \frac 23 + \lambda = 0$$
$$F_y = 2y - \frac 23 + \lambda = 0$$
$$F_z = 2z - \frac 23 + \lambda = 0$$
$$F_{\lambda} = x + y + z -1 = 0$$
Now we get three relation:
$$x = y = z = \frac{\frac {2}{3} - \lambda}{2} = \frac{2-3\lambda}{6}$$
Plugging this into the fourth partial derivative relation we have:
$$x+y+z = 1$$
$$\frac{2-3\lambda}{6} + \frac{2-3\lambda}{6} + \frac{2-3\lambda}{6} = 1$$
$$3(2-3\lambda) = 6$$
$$6-3\lambda = 6$$
$$-3\lambda = 0 \implies \lambda = 0$$
This would imply that:
$$x = y = z = \frac 13$$
But this is the minima, because we have:
$$f(\frac 13, \frac 13, \frac 13) = 0$$
So because this is the only stationary point, the maxima must be on the boundary. And indeed the maxima occurs at these 3 point $$(x,y,z) = (1,0,0), (0,1,0), (0,0,1)$$
$$f(1,0,0) = f(0,1,0) = f(0,0,1) = \frac 23$$
Note that because all these terms are from the same form and we don't have terms of the for $xy$ in the expansion it means that the partial derivatives will have same forms, implying that at the stationary point we'll find the variables will all be the same.
When you check all statinary point in a certain region and you don't get what you want like in this case, then check the boundaries. In this case the region is bounded between $0$ and $1$ so we set one variable at a time to this boundaries and this becomes a 2 variable function which is much easier to maximize.