Perturbation method for a phase coexistence system

Question

The differential equation with $k$ dimensions has the form $\frac{d\mathbf{F}}{d\mathbf{x}} = \mathbf{F}(\mathbf{x},\epsilon)$, where $\mathbf{x}$ is the variable vector and $\epsilon$ is a parameter.  I focus on the steady state, which yields the solution of $\mathbf{F}(\mathbf{x},\epsilon)=0$. Assume that the equation has a perturbation solution $\mathbf{x} = \mathbf{a}\epsilon +\mathcal{O}(\epsilon^2)$ for the parameter $\epsilon\rightarrow 0$ and another perturbation solution $\mathbf{x} = \mathbf{b}_0-\mathbf{b}_1\epsilon^{-1} +\mathcal{O}(\epsilon^{-2})$ for the parameter $\epsilon\rightarrow \infty$.

Now I already know (from numerical simulations) that the above two kinds of solutions i.e., $x_i\rightarrow 0$ for $i\in I$ and $x_j\rightarrow b_0$ for $j\in J$, coexists for the parameter $\epsilon$ of a specific value interval. ($K=\{1,2,..., k\}$ and the sets $I,J\subseteq K$, $I\cap J=\emptyset$).
My question is how to identify the solution sets $I$ and $J$ by perturbation method? or I should consider its Hamiltonian?

Answer 1

Make an Laurent series ansatz $x=\sum_{k\in Z} z_k \varepsilon^{k}$, insert it into the differential equation, expand into a Laurent series and compare coeffiicients to get differential equations for the $z_k$. Then you can analyze a truncated version of this system to find out which coefficients can be set to zero without impairing solvability.