A non-linear wave equation can be written as, \begin{equation} -\phi_{t,t,}(x,t)+\Delta\phi(x,t)=\phi + \sum_{k=2}^\infty g_k \phi^k, \tag{1} \end{equation} the $\phi$ can be exppanded as \begin{equation} \phi=\sum_{k=1}^\infty\varepsilon^k \phi_n \end{equation}
In order to obtain non-trivial solutions of (1) their characteristic scale must also become $\varepsilon$-dependent.What does it mean?
First of all, you want to put the factor of $\varepsilon$ in front of the nonlinear term. Then you may assume that the solution looks like
$$\phi = \phi_0 + \sum_{k=1}^{\infty} \varepsilon^k \phi_k$$
where
$$-\frac{\partial^2 \phi_0}{\partial t^2} + \nabla^2 \phi_0 = \phi_0$$
Now solve to $O(\varepsilon)$, i.e., $\phi = \phi_0+ \varepsilon \phi_1$:
$$-\frac{\partial^2 \phi_0}{\partial t^2} + \nabla^2 \phi_0 + \varepsilon \left ( -\frac{\partial^2 \phi_1}{\partial t^2} + \nabla^2 \phi_1 \right ) = \phi_0+ \varepsilon \phi_1 + \varepsilon g_2 (\phi_0+ \varepsilon \phi_1)^2 \approx \phi_0+ \varepsilon \phi_1 + \varepsilon g_2 \phi_0^2$$
The coefficients of $\varepsilon^0$ all cancel, and equating the coefficients of $\varepsilon^1$, we get
$$ -\frac{\partial^2 \phi_1}{\partial t^2} + \nabla^2 \phi_1 = \phi_1 + g_2 \phi_0^2$$
Because you know $\phi_0$ from the linear solution, you now have an inhomogenoeus, but linear equation for the next tern $\phi_1$. You can solve for this, and repeat for higher-order terms recursively.