I have this system of differential equations which model chemical concentrations in a certain reactions: $$\dot{x}=a-x-\frac{2xy}{1+x^2}\qquad \dot{y}=bx\left(1-\frac{y}{1+x^2}\right)$$ for $a,b >0$ and $x(t),y(t)\geq 0$, and I want to find the steady state of it and then classify it.
I began by finding the only possible steady state, which occurs at $\left(\frac{a}{3},\frac{a^2}{9}+1\right)$, the intersection of $y=\frac{(x^2+1)(a-x)}{2x}$ and $y=1+x^2$.
I then calculated the Jacobian of the system to be $$\left[\begin{matrix} \frac{2(x^2-1)y}{(x^2+1)^2}-1 & \frac{-2x}{x^2+1} \\ \frac{b(x^2-1)y}{(x^2+1)^2}+b & \frac{-bx}{x^2+1} \end{matrix}\right]$$
From here I substituted the values of $x$ and $y$ at the calculated steady state into the Jacobian: $$\left[\begin{matrix} \frac{2(\frac{a^2}{9}-1)}{(\frac{a^2}{9}+1)}-1 & \frac{-2a}{3(\frac{a^2}{9}+1)} \\ \frac{b(\frac{a^2}{9}-1)}{(\frac{a^2}{9}+1)}+b & \frac{-ba}{3(\frac{a^2}{9}+1)} \end{matrix}\right]$$
Now to find the stability of the point $\left(\frac{a}{3},\frac{a^2}{9}+1\right)$, I need to find the eigenvalues of this matrix. I began by trying to calculate this by brute force and it is very messy. Is there an easier way? Have I missed some simplification? Have I made a mistake that is making this a much more complicated calculation than it should be?
I found the characteristic polynomial to be: $$\frac{8a^6b^2}{243\left(\frac{a^8}{2187} + \frac{4a^6}{243} + \frac{2a^4}{9} + \frac{4a^2}{3} + 3\right)} - \frac{8a^4b^2}{27\left(\frac{a^8}{2187} + \frac{4a^6}{243} + \frac{2a^4}{9} + \frac{4a^2}{3} + 3\right)} - \frac{4a^4b^2}{27\left(\frac{a^6}{243} + \frac{a^4}{9} + a^2 + 3\right)} - \frac{8a^3b\lambda}{27\left(\frac{a^6}{729} + \frac{a^4}{27} + \frac{a^2}{3} + 1\right)} + \frac{8a^5b\lambda}{243\left(\frac{a^6}{729} + \frac{a^4}{27} + \frac{a^2}{3} + 1\right)} - \frac{4a^4b^2\lambda}{27\left(\frac{a^6}{243} + \frac{a^4}{9} + a^2 + 3\right)} - \frac{4a^3b\lambda^2}{27\left(\frac{a^4}{81} + \frac{2*a^2}{9} + 1\right)} - \frac{4a^3b\lambda}{27\left(\frac{a^4}{81} + \frac{2*a^2}{9} + 1\right)}$$ and I have no idea how I would solve this for the eigenvalues $\lambda$. There must be a simpler way!
i am going to set $a = 3$ see what happens. i will see if by scaling time if this can be done. i will make a change of variable $$x = 1 + u, y = 2 + v.$$ then
$$u' = 3 - (1+u) - \frac{2(1+u)(2 + v)}{1 + (1 + u)^2 } = \frac{(2-u)(2+2u+u^2)-2(2+2u+v+uv)}{2+2u+u^2} = \frac{(4-2u+4u+\cdots) -(4+4u+2v+\cdots)}{2+2u+u^2} = -u -v+\cdots$$ and
$$ v' = b(1+u)\left(1 - \frac{2+v}{2+2u+u^2}\right) = b(1+u)\frac{2u-v+u^2}{2+2u+u^2}= b(2u-v+\cdots)$$ the linear stability is determined by $$u' = -u - v, \, v' = 2bu - bv. $$ the characteristic equation is $$\lambda^2 + (b+1) \lambda + 3b = 0, \lambda = \frac{-b-1 \pm \sqrt{b^2-10b + 1}}2 $$ so if $b > 5 + 2\sqrt 6,$ then the real part of $\lambda < 0$ implying a stable spiral. if $\lambda < 5 + 2\sqrt 6,$ then because roots are negative, we have a stable focus.