I have the following two equations with complex $\alpha, \beta$: $$0 = -i \omega_m \beta - ig|\alpha|^2 - \frac{\gamma}{2}\beta$$ $$0 = \alpha(i\Delta - \frac{\kappa}{2}) - i\epsilon - ig\alpha(\beta + \beta^*)$$ The other parameters $g, \omega_m, \Delta, \kappa, \gamma, \epsilon$ are real.
To solve this system for $\alpha$, I solved the first equation for $\beta$ and got $$\beta = \frac{g|\alpha|^2}{\frac{i\gamma}{2} - \omega_m}$$
I then inserted this into the second line and obtained $$\alpha\left(i\Delta - \frac{\kappa}{2} + \frac{2ig^2\omega_m|\alpha|^2}{\frac{\gamma^2}{4} + \omega_m}\right) = i\epsilon$$However, in this equation, both $\alpha$ and $|\alpha|^2$ appear, and I don't know how to solve this.
The last equation is of the form $\,\alpha\left(u + i\left(v + w\, |\alpha|^2\right)\right) = i \epsilon\,$ with $\,u,v,w,\epsilon \in \mathbb R\,$.
Taking squared magnitudes on both sides:
$$ |\alpha|^2\left(u^2 + \left(v+ w\,|\alpha|^2\right)^2\right) = \epsilon^2 $$
This is a cubic equation with real coefficients in $\,|\alpha|^2\,$ which can be solved either algebraically or numerically. Substituting the positive real solution(s) $\,|\alpha|^2\,$ in the first equation of the system gives $\,\beta\,$, then substituting the respective $\,\beta\,$ in the second equation gives $\,\alpha\,$.