Let us look at an example, the boundary layer problem $$\epsilon y'' + (1+x) y'+ y = 0$, \quad y(0) = y(1) = 1.$$ We will substitute $y(x) = \exp\left[\frac{1}{\delta} \sum_{n=0}^\infty \delta^n S_n(x) \right]$ into the equation and divide by $y$, and we get $$\left[\frac{\epsilon}{\delta} S_0'' + \epsilon S_1 '' + \cdots + \frac{\epsilon}{\delta^2} S_0'^2 + 2\frac{\epsilon}{\delta}S_0' S_1' \cdots\right] + \left[\frac{1}{\delta}(1+x) S_0' + \cdots\right] + 1 = 0 \quad (\star)$$
Two largest terms $\frac{\epsilon}{\delta^2} S_0'^2$ and $\frac{1}{\delta}(1+x) S_0'$ must be balanced, so we take $\delta = \epsilon$ and at $O(1/\epsilon)$, we have $$S_0'^2 + (1+x)S_0' \sim 0.$$
There are two cases, $S_0' \sim 0$ or $S_0' \sim -(1+x)$.
For the first case, we get $S_0 = A$ a constant, and look at $(\star)$ at $O(1)$ level, we get $$\frac{1}{\epsilon} S_0'^2 + \frac{1}{\epsilon} (1+x) S_0' + S_0'' + 2S_0' S_1' + (1+x) S_1' +1 \sim 0$$ which reduces to $$(1+x)S_1' \sim -1$$ so $S_1 \sim -\log(1+x) + B$. Now this will give a solution $$y_1 = \exp\left[\frac{1}{\epsilon}S_0 + S_1 + O (\epsilon)\right] \sim \exp\left[\frac{1}{\epsilon}A + \Big(-\log(1+x) + B\Big)\right] \text{ as }\epsilon \rightarrow 0.$$
The next step, the book says the above is $\sim \frac{C}{1+x}$ where $e^{\frac{1}{\epsilon}A}$ is obsorbed into the constant $C$, but how is this possible since this term is not small as $\epsilon$ goes to zero.
For completeness, I will add the calculation for the second case where $S_0' = -(1+x)$. Then we have $S_0 = -x^2/2 - x + B$. Again looking at $O(1)$ terms, $$\frac{1}{\epsilon} S_0'^2 + \frac{1}{\epsilon} (1+x) S_0' + S_0'' + 2S_0' S_1' + (1+x) S_1' +1 \sim 0$$ this reduces to $$-1 -2(1+x)S_1' + 1 \sim 0$$ so $S_1' \sim 0$ and $S_1 = C_1$. So we get a second solution $$y_2 \sim \exp\left[\frac{1}{\epsilon}\left( -\frac{x^2}{2} - x + B \right) + C_1 \right]$$
Again, the book claims the above is $\sim > D\exp\left[\frac{1}{\epsilon}\left( -\frac{x^2}{2} - x \right) > \right]$, so the constant $e^{B/\epsilon}$ is also absorbed into the constant $D$.
You do the integration for fixed $ϵ$, or in other words, the integration constants $A,B,C$ etc. can depend on $ϵ$, they are only constant relative to $x$. The constants are fixed by the boundary conditions.
Relative to https://math.stackexchange.com/a/2807460/115115 where $f(x)=x+1$ and the anti-derivatives with $F(0)=0=G(0)$ are $F(x)=\int f dx=\frac12x^2+x$, $G(x)=\int \frac1fdx=\ln(x+1)$, the solution form is $$ y(x)\approx Ae^{-G(x)}+\frac{B}{f(x)}e^{-\frac1εF(x)+G(x)}=\frac{A}{x+1}+Be^{-\frac1{2ε}(x^2+2x)}. $$ To satisfy the initial conditions \begin{align} 1=y(0)&=A+B\\ 1=y(1)&=\frac A2+Be^{-\frac3{2ε}} \end{align} we get approximatively $A=2$ and $B=-1$, as $ e^{-\frac3{2ε}}\in O(ε^k)$ for any $k>0$. The solution that closely approximates the solution and the boundary conditions is thus $$ y(x)\approx \frac{2}{x+1}-e^{-\frac1{2ε}(x^2+2x)} $$