R.<u,v,x,y,z> = PolynomialRing(QQ, order='lex')

p = (x-u)^2 + (y-v)^2 + z^2 - 1/10
dpdt = -v*diff(p, u) + u*diff(p, v)

time = u^2 + v^2 - 1

I = ideal(p, dpdt, time)
G = I.groebner_basis()
eqn = G[-1]

print(eqn)
# x^4 + 2*x^2*y^2 + 2*x^2*z^2 - 11/5*x^2 + y^4 + 2*y^2*z^2 - 11/5*y^2 + z^4 + 9/5*z^2 + 81/100


# plotting

var('x y')
eqn = x^4 + 2*x^2*y^2 + 2*x^2*z^2 - 11/5*x^2 + y^4 + 2*y^2*z^2 - 11/5*y^2 + z^4 + 9/5*z^2 + 81/100
implicit_plot3d(eqn, (-2, 2), (-2, 2), (-2, 2))