Consider the non-linear system $$ \frac{d\mathbf{x}}{dt} =f\left(\mathbf{x},\mathbf{u} \right) $$ where $f: \mathbb{R}^{n+2}\mapsto\mathbb{R}^n$ and parameters $\mathbf{u}\in\mathbb{R}^2$. I am interested in performing a bifurcation analysis around a bifurcation point $\mathbf{u}=\mathbf{u}_c$. Assume the homogeneous steady-state is hyperbolic (that is, all the eigenvalues of the jacobian at a such point have non-zero real parts). How can I perform, for example, a multiple-scale method to analyze the dynamics around this bifurcation?
For example, take the system $$ \frac{dx_i}{dt}=\frac{(y_{i-1}+y_{i+1})^2}{u_1+(y_{i-1}+y_{i+1})^2}-x_i\\ \frac{dy_i}{dt}=\frac{1}{1+u_2x_i^2}-y_i $$ for $1\leq i\leq n$. I can linearize the system around a homogenous steady-state $(x^*,y^*)$ and decouple the system using a Fourier transform, so that, for a small perturbation $(\tilde{x},\tilde{y})$, $$ \frac{d}{dt}\begin{pmatrix}\tilde{x}_i\\\tilde{y}_i\end{pmatrix}\simeq \mathbf{J}(u_1,u_2)\begin{pmatrix}\tilde{x}_i\\\tilde{y}_i\end{pmatrix} $$ where $\mathbf{J}(u_1,u_2)$ is the Jacobian matrix of the original non-linear system plus some decoupling matrix, evaluated at $(x^*,y^*)$. The solution to the linearized system is a combination of exponentials, which parameters depend on the eigenvalues of $\mathbf{J}$. I would like now to understand the bifurcation at the critical values of the bifurcation parameters $(u_1,u_2)$, that is, where the linearized solution changes its behaviour. Any thoughts?