The goal is to find where the normal of the two surfaces are parallel, that will be the maximum and minimum (maybe saddle as well) if I remember my analysis course.
Using lagrange multipliers this is expressed as a number of equations of which has to be fulfilled to reach the goal above.
The equations consist of the partial derivates with respect to x,y,z of F(x,y,z) and the constraint G(x,y,z)=x^2 + y^2 + z^2 -1
To simplify F_x means F differentiated with respect to x, and m is the lagrange multiplier. With that notation you get these equations:
F_x=m*G_x
F_y=m*G_y
F_z=m*G_z
You also have the constraint:
x^2 + y^2 + z^2=1, which also must be fulfilled.
So in total you have 4 equations from which you can find a number of points and the lagrangian multiplier.
In your case you get:
2x*y^2*z^2=2m*x
2y*x^2*z^2=2m*y
2z*y^2*x^2=2m*z
x^2 + y^2 + z^2=1
The three first implies that: y^2=x^2=z^2 => y=x=z
Inserting that for each variable in the last equation
you get:
3x^2=1=>x=(+-)*(1/3)^(1/2)
By using y=x=z you get two points
-((1/3)^(1/2), (1/3)^(1/2), (1/3)^(1/2))
and
((1/3)^(1/2), (1/3)^(1/2), (1/3)^(1/2))
I'm a little unsure of how to check if they are min, max or saddle. For two variables it's pretty straightforward. Maybe you could develop the taylor series around the points and see what happens if you increase/decrease the values of the points slightly.