#Try: S({x+y-2,x+2*y-3},{x,y});

S:=proc(eq,var) local v,i,a,eq1,var1,S1:
if eq={} then  RETURN({seq(v=v, v in var)}): fi:
if var={} and eq<>{} then  RETURN(FAIL): fi:
v:=var[1]:

for i from 1 to nops(eq) while coeff(eq[i],v,1)=0 do od:

if i=nops(eq)+1 then S1:=S(eq,var minus {v}):  

if S1=FAIL then   RETURN(FAIL):  else RETURN(S1 union {v=v} ):  fi: fi:

a:=-coeff(eq[i],v,0)/coeff(eq[i],v,1):

eq1:=expand(subs(v=a,eq)) minus {0}: var1:=var minus {v}:

S1:=S(eq1,var1): if S1=FAIL then  RETURN(FAIL): fi:

S1 union {v=subs(S1,a)}: 
end: