Multiplying both sides of the first equation by $z$, of the second by $x$, and of the third by $y$, we get
$$x^2z=z^2y+z, \quad y^2x=x^2z+2x,\quad z^2y=y^2x+4y.$$
Adding up and cancelling, we get
$$2x+4y+z=0.$$
Similarly,
$$x^2y=y^2z+y, \quad y^2z=z^2x+2z, \quad z^2x=x^2y+4x,$$
giving
$$4x+y+2z=0.$$
To finish, use the linear equations to express $y$ and $z$ in terms of $x$, and substitute in $x^2=yz+1$.
Your objective function $c^\top x$ seems wrong. Just use
\begin{matrix}
\min & 0^\top x \\
\text{w.r.t.} & A x = b
\end{matrix}
and your LP solver should return no solution if there is no solution for $Ax = b$ or otherwise one of the solutions.
Example using lpSolve from R:
If not done already, install lpSolve:
> install.packages("lpSolve")
> install.packages("lpSolveAPI")
Load the library:
> library(lpSolveAPI)
Create a new model in $5$ unknowns and store it in the variable lprec:
> lprec<-make.lp(0,5)
Set the objective function t0 $c = 0$:
> set.objfn(lprec, c(0,0,0,0,0))
Now add the constraints, thus the rows of $A x = b$:
> add.constraint(lprec, c(1,1,0,1,0), "=", 10)
> add.constraint(lprec, c(0,1,1,1,0), "=", 12)
> add.constraint(lprec, c(1,1,1,0,0), "=", 9)
> add.constraint(lprec, c(1,0,1,1,0), "=", 11)
> add.constraint(lprec, c(1,1,0,0,1), "=", 15)
> add.constraint(lprec, c(1,1,0,0,0), "=", 5)
> add.constraint(lprec, c(0,1,1,0,0), "=", 7)
Have a look at the model:
> lprec
Model name:
C1 C2 C3 C4 C5
Minimize 0 0 0 0 0
R1 1 1 0 1 0 = 10
R2 0 1 1 1 0 = 12
R3 1 1 1 0 0 = 9
R4 1 0 1 1 0 = 11
R5 1 1 0 0 1 = 15
R6 1 1 0 0 0 = 5
R7 0 1 1 0 0 = 7
Kind Std Std Std Std Std
Type Real Real Real Real Real
Upper Inf Inf Inf Inf Inf
Lower 0 0 0 0 0
Now solve:
> solve(lprec)
[1] 0
And check the found optimal value:
> get.objective(lprec)
[1] 0
Then check for the found feasible solution
> get.variables(lprec)
[1] 2 3 4 5 10
Best Answer
Let x+y=u and xy=v
Second equation becomes u(v-1)=3
First equation is ${(xy)}^2+(x+y)^2+1-2xy=10$ which is $v^2+u^2+1-2v=10 \Rightarrow (v-1)^2+u^2=10$
Eliminate v-1 from the two equations to get $\frac{9}{u^2}+u^2=10$
Solve the quadratic to get u=$\pm1$ and $\pm3$ then solve for x and y