We plot Steiner's surface, given implicitly by the equation: $$x^2y^2+y^2z^2+z^2x^2=xyz$$ First, we initialize the variables and define the function

def f(x,y,z):
    return x^2*y^2+y^2*z^2+z^2*x^2-x*y*z
x,y,z=var('x y z')

Now we plot it:

steiner = implicit_plot3d(f(x,y,z)==0,(x,-.5,.5),(y,-.5,.5),(z,-.5,.5),plot_points=100,opacity=0.75)
show(steiner,viewer='tachyon')

It's instructive to find the critical points, so let's compute the derivative of $f$

df = [diff(f(x,y,z),w) for w in [x,y,z]]
df

[2*x*y^2 + 2*x*z^2 - y*z, 2*x^2*y + 2*y*z^2 - x*z, 2*x^2*z + 2*y^2*z - x*y]

Let's solve for when this vanishes:

CritPoints = solve(df,x,y,z,solution_dict=True)
CritPoints

[{x: 0, z: r4, y: 0},
 {x: r5, z: 0, y: 0},
 {x: 0, z: 0, y: r6},
 {x: 1/4, z: -1/4, y: -1/4},
 {x: -1/4, z: 1/4, y: -1/4},
 {x: 1/4, z: 1/4, y: 1/4},
 {x: -1/4, z: -1/4, y: 1/4}]

Notice that first three entries are the $x,y,$ and $z$ axes. This explains why something funny seems to happen in the plot of the surface along those axes.

Ok. Now, let's use this list of critical points to find the critcal values of $f$:

for p in CritPoints:
    print f(x,y,z).subs(x=p[x],y=p[y],z=p[z])

0
0
0
-1/256
-1/256
-1/256
-1/256

So the set of regular values are $\mathbb{R}\setminus\{0,-1/256\}$ Let's try plotting the surface again, but now at 1/512

aSteiner = implicit_plot3d(f(x,y,z)==1/512,(x,-.5,.5),(y,-.5,.5),(z,-.5,.5),opacity=0.75)
show(aSteiner,viewer='tachyon')

This is indeed a manifold now, but it is not Steiner's surface.