Plot intersections of two relations

edit

Since you have two equations and three variables, first solve the system in terms of one variable:

var('x y z')
f=solve( [ x^2+(y-1)^2+(z-1)^2==1, x^2+y^2+z^2==1 ], y, z )

This gives solutions as a list:

[[y == -1/2*sqrt(-2*x^2 + 1) + 1/2, z == 1/2*sqrt(-2*x^2 + 1) + 1/2], [y == 1/2*sqrt(-2*x^2 + 1) + 1/2, z == -1/2*sqrt(-2*x^2 + 1) + 1/2]]

You can then either copy the solutions into a 3d parametric plot,

p = parametric_plot3d( [ x, -1/2*sqrt(-2*x^2 + 1) + 1/2, 1/2*sqrt(-2*x^2 + 1) + 1/2 ], (x,-1,1) )
p += parametric_plot3d( [ x, 1/2*sqrt(-2*x^2 + 1) + 1/2, -1/2*sqrt(-2*x^2 + 1) + 1/2 ], (x,-1,1) )

or substitute the solutions through the list indices:

p = parametric_plot3d( [ x, y.subs(f[0][0]), z.subs(f[0][1]) ], (x,-1,1) )
p += parametric_plot3d( [ x, y.subs(f[1][0]), z.subs(f[1][1]) ], (x,-1,1) )

and then show the combined plot. Here's an example including a unit sphere in cyan and the parametric solution in blue, with a slightly different output method.