%% multipol example
% create the neccessary number of variables for your problem (2 here)
x = create_vars(2)
% write the equations
eqs(1) = 25*x(1)*x(2) - 20*x(2) - 15*x(1) + 12
eqs(2) = x(2)^2 + x(1)^2 -1
%% obtain the matrix of coefficients and vector of monomials (choose ordering)
[c mon] = polynomials2matrix(eqs,'plex')
%% add new equations, here chosen randomly
eqs(3) = x(1)*eqs(1)
eqs(4) = x(2)*eqs(2)
%% obtain new coefficient matrix and vector of monomials
[c mon] = polynomials2matrix(eqs,'plex')
%% check if GJ elimination gives you an equation in single variable
cr = rref(c)
% the last row corresponds to y^3-0.6y^2-0.36y+0.216
%% find the solutions for y
yr = real(roots(cr(end,end-3:end)))
% now you would need to backsubstitute and find corresponding x's and then
% verify all pairs of (x,y) solutions to see if they satisfy the original
% equations