I have a system of coupled differential equations (an example )in the form of, x &#x203

I have a system of coupled differential equations (an example )in the form of,
$x^{″}+a{x}^{\prime }+bx-<!--\; -\; -->cy=0$
$y^{″}-<!--\; -\; -->a{y}^{\prime }+by-<!--\; -\; -->cx=0$
The solution to the above system looks like,
$x=A{e}^{w_{1}t}+B{e}^{w_{2}t}+C{e}^{w_{3}t}+D{e}^{w_{4}t}$
$y=E{e}^{w_{1}t}+F{e}^{w_{2}t}+G{e}^{w_{3}t}+H{e}^{w_{4}t}$
the frequencies $w_1$,$w_2$,$w_3$ and $w_4$ are functions of a,b and c My goal is to determine that values of a,b and c such that all w are Real. However, my original system of differential equations has nonlinear terms. I managed to derive the solution but in power expansion form or,
$x=\underset{i=0}{\overset{n}{\sum <!--\; \sum \; -->}}\frac{a_i}{i!}{t}^i$
$y=\underset{i=0}{\overset{n}{\sum <!--\; \sum \; -->}}\frac{a_i}{i!}{t}^i$
I need to determine the frequencies so that I can find the values of a,b and c such that all w are Real. Is it possible to derive w from limited power expansion (say 8th term) of x and y?
My first attempt:
I did a test of my method by taking Fourier Sine Transform of analytical solution which give me the answer in form of $\frac{f(w)}{g(w)}$ where w is the frequency. Then solved the g(w) for w which gives me four answers for w which are the frequencies $w_1$,$w_2$,$w_3$ and $w_4$. However when I attempted with Fourier Sine Transform of series solutions, the answer is different. This is due to the limited power expansion. Is here a way to improve this?
My second attempt:
I linearize the nonlinear terms of differential equations and used the matrix to calculate Determinant of (A-wI) where A is the matrix of the system of differential equation. I managed to calculate the values of a,b and c but they were incorrect because of linear terms I made.
A last method that I am considering is Monte Carlo. Is it possible to get a frequency equation from limited power expansion of differential equation solution? Any ideas for other methods that I missed?