restart
n=2 --specify n here

QQ[x_1..x_n]
scan(n, i -> u_(i+1) = random(QQ)) --specify u here
I = ideal((x_1^2 + x_2^2 + x_2)^2 - (x_1^2 + x_2^2)) 
--specify generator of I(X) here

c = codim I
J = transpose jacobian I
xMinusU = matrix {apply(n, i -> x_(i+1)-u_(i+1))}
AJ = xMinusU || J
CritCond = I + minors(c+1,AJ)

Sing = radical ideal singularLocus I
Crit = saturate(CritCond,Sing)

dim Crit, degree Crit