How to eliminate secular terms for perturbed non-oscillatory equations? Even in a linear second order equation like $x''+x'+\epsilon x=0$ the standard asymptotic expansion has a secular term already in the first order of $\epsilon$, namely 
$$x(t)=a_0+b_0e^{-t}+\epsilon(a_1+a_0t+b_1e^{-t}-b_0te^{-t})+O(\epsilon^2). 
$$
But the exact solution is obviously bounded uniformly in $t,\epsilon\geq0$, the expansion therefore is not valid for large times. The usual approach to eliminating secular terms is to use multiple time scales, but all perturbation theory texts I looked at (Holmes, Hunter, Kevorkian-Cole, Verhulst) only consider cases where the unperturbed equation is oscillatory. The Poincare-Lindstedt method or averaging that are used only make sense when there are oscillations. But the above equation is obviously non-oscillatory for small $\epsilon$. I am interested in non-linear equations that behave similarly to the linear example above, e.g. 
$$x''+x'\Big(1-\frac32 \frac{x'}{x}\Big)+\epsilon x^3=0.$$ They come up when approximately realizing mechanical constraints by viscous friction, which explains the absence of oscillations. 

The issue seems to be that bounded solutions that are not exponentially stable produce secular terms even without oscillations, so one has all the pain of secular terms with no benefit of averaging. For instance, I suspect that for small $\epsilon$ solutions to the above nonlinear equation have a limit when $t\to\infty$, which depends on the initial values, it would be nice to approximate its dependence on them and $\epsilon$. Is there something of this sort done in the literature? References are appreciated. 

 A: This is probably more a comment than an answer, but too long for the former.
A different approach might be to use normal hyperbolicity theory to first find a series expansion of an invariant manifold for your system. Then, you might be able to express solutions as a combination of an outer solution converging to the invariant manifold, and an inner/boundary layer solution on the invariant manifold (I'm not sure my use of terminology is completely correct here).
In your case, the invariant manifold at $\epsilon = 0$ is $M = \{ \dot{x} = 0 \}$. This is exponentially attractive and the dynamics on $M$ is trivial. This implies that $M$ is normally hyperbolic and then allows you to expand the persistent $M_\epsilon$ in $\epsilon$ explicitly, order by order. Then, you can express the dynamics on $M_\epsilon$ using this expansion.
I replied a little later than I first saw this question, because I was just finishing an article that shows this idea, see http://arxiv.org/abs/1603.00369. I'm also studying systems with friction and show that in a limit of sending friction to infinity one obtains nonholonomic dynamics, but also one can do a (singular, not regular) expansion of the dynamics. The Chaplygin sleigh example (including section 7.1) might be instructive.
Added notes on non-compactness
First of all, even on a compact space, the solution curves of the approximate system may not converge uniformly in time to solutions of the nonholonomic system. As counterexample see e.g. the nonholonomic constraint realization example in section 3 of my preprint, but simply compactify the space away from the original circle: any approximate solution will still drift a given finite distance away from the circle, if you wait long enough. However, the solution to the approximate system is at all times an $\mathcal{O}(\epsilon)$-pseudo solution of the nonholonomic system, that is, on any finite time interval such solutions converge uniformly to a solution of the nonholonomic system. I don't know of a geometric condition for when this would improve to uniform convergence for all positive time, i.e. when the drift would not accumulate. Since this drift is essentially given by the reaction force along solutions of the nonholonomic system, you'd need to know these.
For non-compact spaces another reason that the solutions of the approximate system do not converge uniformly may be because of non-uniform exponential stability of the invariant manifold. That may lead to the perturbed invariant manifold to not be uniformly $\mathcal{O}(\epsilon)$-close to the unperturbed one. This you can check by inspecting normal attraction rate and uniform $C^1$-boundedness of the vector field near the invariant manifold. In your example the attraction rate is
$$
-\Bigl(1-\frac{3}{2}\frac{x'}{x}\Bigr),
$$
which is uniformly $1$ on $M = \{ x' = 0 \}$, but it blows up at $x = 0$, when $x' \neq 0$ but small. Hence, for solution curves passing by that point you may not even get convergence in the $\mathcal{O}(\epsilon)$-pseudo solution sense.
