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.