Consider the following differential equation
$$ y''=-y + \alpha y |y|^2, $$
where $y=y(x)$ is complex in general and $\alpha$ is a real constant such that the second term is small compared to $y$ ($||^2$ is the absolute square). Numerically I find that a good approximation for this is found by taking the solution for $\alpha=0$, which is $ae^{ix}+be^{-ix}$ (with $a$ and $b$ determined by the initial conditions) and applying a little shift to the exponents: $ae^{i(1+\Delta_a)x}+be^{-i(1+\Delta_b)x}$, where $\Delta_a=-\alpha(|a|^2+2|b|^2)/2$ and $\Delta_b=-\alpha(2|a|^2+|b|^2)/2$.
I want to find a simple justification for this approximate solution. My attempt so far is this: assuming the solution $ae^{i(1+\Delta_a)x}+be^{-i(1+\Delta_b)x}$, and inserting it in the differential equation, we find (neglecting the $\Delta^2$ terms, and expanding the absolute square):
$$\Delta_a ae^{i(1+\Delta_a)x}+\Delta_b be^{-i(1+\Delta_b)x}=-\alpha/2 \bigg [ (|a|^2+2|b|^2)ae^{i(1+\Delta_a)x} + (2|a|^2+|b|^2)be^{-i(1+\Delta_b)x} +ab^*ae^{i(3+2\Delta_a+\Delta_b)x}+ba^*be^{-i(3+\Delta_a+2\Delta_b)x} \bigg ]$$
By matching the exponentials of same power, we directly find the expression of the $\Delta$s. However there are those two (roughly-) third power terms left, that seem to indicate at first that the assumed form is not correct (because these are not small compared to the first two terms). And yet the $\Delta$s found in this way give a very good approximation. (You can choose initial conditions such that the last two terms are small compared to the first two, but my approximate solution is good even when this is not the case.)
How can we somehow disregard the last two terms? Or can we justify the $\Delta$s in a different way?
Yes, due to the non-linear perturbation the frequencies will mix, and the base frequency also gets perturbed. You get terms with frequency about 3 on the right side that have coefficients of order $O(\alpha^1)$. Neglecting them would mean to reduce the approximation order to zero and you would be back at the unperturbed equation. Thus you have to compensate these terms on the right side by adding corresponding frequency 3 terms to $y$ and thus to the left side.
This will also lead to terms of up to frequency 7 on the right side, but the coefficients there will have order 2 or 3 in $\alpha$, thus they can be ignored if the whole equation is considered only to first order, with a residual that is $O(\alpha^2)$.
So $$ y = a_3e^{i(3+2Δ_a+Δ_b)x}+a_1e^{i(1+Δ_a)x} +b_1e^{-i(1+Δ_b)x}+b_3e^{-i(3+Δ_a+2Δ_b)x} $$ and for the terms it makes more sense to use a table from the start \begin{array}{c|c|c|c} \text{frequeny}&\text{combination} (y,y|\bar y)&\text{coeff. left}&\text{coeff. right}\\ \hline 1+Δ_a&(1,1|1)&-2Δ_a&|a_1|^2 a_1 \\ &2(1,-1,-1)&&2|b_1|^2a_1\\ &2(3,-1|1)&&2a_3b_1\bar a_1\\ \hline -1-Δ_b&(-1,-1|-1)&-2Δ_b&|b_1|^2 b_1\\ &2(1,-1|1)&&2|a_1|^2b_1\\ &2(-3,1|-1)&&2b_3a_1\bar b_1 \\ \hline 3+2Δ_a+Δ_b&(1,1|-1)&-8-3(2Δ_a+Δ_b)&a_1^2\bar b_1 \\ &2(3,1|1)&&2a_3|a_1|^2 \\ &2(3,-1|-1)&&2a_3|b_1|^2 \\ \hline -3-Δ_a-2Δ_b&(1,1|-1)&-8-3(Δ_a+2Δ_b)&b_1^2\bar a_1 \\ &2(-3,1|1)&&2b_3|a_1|^2 \\ &2(3,-1|-1)&&2b_3|b_1|^2 \\ \end{array} There are further terms in these frequencies that combine two or more factors of frequency 3, which surpasses the approximation order in the coefficients.
So in comparing coefficients of the same frequency we get \begin{align} -2Δ_aa_1 &= \alpha ((|a_1|^2+2|b_1|^2) a_1 + 2a_3b_1\bar a_1)\\ -2Δ_bb_1 &= \alpha ((2|a_1|^2+|b_1|^2) b_1 + 2b_3a_1\bar b_1)\\ -(8+3(2Δ_a+Δ_b))a_3 &=\alpha (a_1^2\bar b_1 + 2a_3(|a_1|^2 + |b_1|^2))\\ -(8+3(2Δ_a+Δ_b))b_3 &=\alpha (b_1^2\bar a_1 + 2b_3(|a_1|^2 + |b_1|^2))\\ \end{align} So in the lowest order, \begin{align} a_3 = -\frac{\alpha}{8}a_1^2\bar b_1\\ b_3 = -\frac{\alpha}{8}b_1^2\bar a_1\\ \end{align} which then reduces the first two equations in only their first-order terms to the form you already found.
Using all the terms in the table would give you the approximation in $O(\alpha^3)$, but would need the frequency 5 terms to balance. In conclusion $$ y = -\frac{\alpha}{8}a_1^2\bar b_1e^{i(3+2Δ_a+Δ_b)x}+a_1e^{i(1+Δ_a)x} +b_1e^{-i(1+Δ_b)x}-\frac{\alpha}{8}b_1^2\bar a_1e^{-i(3+Δ_a+2Δ_b)x} $$ where $$ Δ_a=-\frac{\alpha}{2}(|a_1|^2+2|b_1|^2)\text{ and } Δ_b=-\frac{\alpha}{2}(2|a_1|^2+|b_1|^2). $$