Define ZeroDimRad(I)
If Type(I)<>IDEAL Then Return Error("Input must be an Ideal"); EndIf;
Ind:=Indets();

Sqf:=[];
  Foreach X In Ind Do
	Elimid:=Elim(Diff(Ind,[X]),I);
	If Elimid=Ideal(0) Then Return Error("Not a Zero Dimensional Ideal");EndIf;
	Gen:=Gens(Elimid);
	G:=Gen[1]/LC(Gen[1]);
	Append(Sqf,SqFree(G));
  EndForeach;
  J:= I + Ideal(Sqf);

Return Minimalized(J);

EndDefine;

 OUTPUTS
a) Use S::=Q[x,y];
   I1:=Ideal(x^3, x^2y + x, y^2); 
   ZeroDimRad(I1);
   Ideal(x, y)
-------------------------------
b) Use Z3::=Z/(3)[x,y];
   I2:=Ideal(x^27+y^27+1,x^18-x^9y^9-x^9); 
   ZeroDimRad(I2);
   Ideal(x, y + 1)
-------------------------------

c) Use R::=Q[x,y,z]; 
   I3:=Ideal(y^4+2y^3+y^2, xy^2+y^3+xy+y^2,y^2z+yz,z^2,xz+yz,x^2+2xy+y^2);
   ZeroDimRad(I3);
   Ideal(x^2 + 2xy + y^2, x^2 - x, y^2 + y, z)
-------------------------------

d)  Use Z7::=Z/(7)[x,y,z];
    I4:=Ideal(y^98+y^49, x^49+y^49,z^7);
    ZeroDimRad(I4);
    Ideal(x^49 + y^49, x^2 - x, y^2 + y, z)
-------------------------------