I have a system of equations of the form $$ \begin{aligned} a &=z+\frac{x}{y^2} &(1)\\ b &=x-2c\frac{x}{y}+ac^2 &(2)\\ c &=y-\frac{z}{x}y^3 &(3), \end{aligned} $$ where $a>0$, $b\geq 0$, and $c>0$. The solution $(x,y,z)$ represents parameters of a physical system that must be positive to make sense (so nonpositive or complex solutions make no sense for my problem).
To solve the system, I first solved $(2)$ for $x$ to get $$ x=-\frac{y \left(b-a c^2\right)}{2 c-y}. $$ Substituting this solution for $x$ into $(1)$ and solving the result for $z$ gives $$ z=\frac{-a c^2+2 a c y-a y^2+b}{y (2 c-y)}. $$ Substituting these forms of $x$ and $z$ into $(3)$ then gives (after some basic algebraic simplification) $$ \tag{*} c = y + ((b - a (c - y)^2) y)/(b - a c^2). $$ So $(*)$ is a cubic equation in $y$, which has three potential solutions. Once I know which of these solutions is the one I need, I instantly obtain the final solutions for $x$ and $z$ as they are in terms of $y$. But how do I know which of these solutions is the one I am looking for?
What I have tried so far:
I tried having Mathematica solve $(*)$ for me, which gave me three solutions. Call them $y_1$, $y_2$, and $y_3$. Substituting in some possible values for $a$, $b$, and $c$, I found only the $y_1$ solution was positive while the other two were complex. Does this mean that $y_1$ is the solution I am looking for? My concern is maybe for different $a$, $b$, and $c$ it may be possible that $y_1$ ends up being nonpositive or complex and instead $y_2$ or $y_3$ gives the only positive solution.
To better understand the nature of the roots I computed the polynomial discriminant of $(*)$, which after some simplifications gives $$ \Delta=\frac{a \left(3 a b c^2+3 \left(b-a c^2\right)^2+29 b^2\right)}{\left(b-a c^2\right)^3}. $$ The numerator is never negative so I concluded $\operatorname{sign}(\Delta)=\operatorname{sign}(b-a c^2)$; however, I am unable to make any conclusions about $\operatorname{sign}(b-a c^2)$ as it could be positive or negative. If it were always negative then I would know $y_1$ is the solution I'm looking for as such a result would guarantee(?) that $y_2$ and $y_3$ are a conjugate pair.
So where I'm at right now is for all the values of $a$, $b$, and $c$ that I've tried, $y_1$ is the only solution that makes sense; however, I am not sure if that implies $y_1$ is the correct solution for all possible $a$, $b$, and $c$. Thoughts?
This is quite a tricky problem. The comment by Robert Israel points us in the right direction:
When we plot (3) for $z,x$ constant, we see the following:
(Ignore the label on the vertical axis).
Where the vertical lines are different values for $c$ and the curved plot is $y-\frac{z}{x}y^3$. You see that the solutions of this equation are related to the solutions of the equation $$y-\frac{z}{x}y^3=0,$$ which can be solved to give $y=0$ or $y=\pm\sqrt{\frac xz}$. So we have two solutions $y>0$ and one solution $y<0$. Let these three solutions be $y_1<y_2\le y_3$.
In particular, if $c$ is too big, there is only one negative solution and two complex solutions.
We can thus conclude that If there are two complex solutions, then there are no positive solutions.
This already is an answer to your question: You should not take your $y_1$, but instead one of the other two solutions. This is what Robert Israel meant.
Furthermore, if $\Delta<0$, there is only one solution for $y$, and from this plot we see that it will be $<0$. So we require $\Delta\ge0$ and thus $b-ac^2\ge0$. (In particular, there are no solutions of interest to you when $b<ac^2$).
Combining this with (2) and $x>0$ we see that $y>2c$.
Combining this with (3) we get $y>2(y-\frac zxy^3)$ which (because $y>0$) gives $1>2-1\frac zxy^2$ which gives $y>\sqrt\frac x{2z}$.
Furthermore, we have that $y_2\le y_m$ where $y_m$ is the $y$-value where $y-\frac zxy^3$ takes a maximum. By setting the derivative to $0$ we find that $y_m=\sqrt{\frac x{3z}}<\sqrt{\frac x{2z}}$. So in particular, the value $y_2$ does not work either (As it will give $x<0$) so you will need the biggest value of $y$ that solves your equation.