using HomotopyContinuation, LinearAlgebra
@var x1 x2
u1 = 2; u2 = 1;
f = (x1^2 + x2^2 + x2)^2 - (x1^2 + x2^2);
AJ = [[x1 - u1; x2 - u2] differentiate(f, [x1; x2])];
solve([f; det(AJ)])