perturbation method for $y' - (y-1)^2 = \epsilon \frac{y^2}{x^2}$ and $y(1) = 1$. The question follows this exercise in Bender and Orszag.
We would like to find the asymptotic for 
$$y' - (y-1)^2 = \epsilon \frac{y^2}{x^2}\text{ and } y(1) = 1$$
Dropping the epsilon term, we get two solutions 
$$y=1 \text{ or } y = 1- \frac{1}{x-c}$$
and the answers there said we will "add perturbation terms to either $y$ or $(y-1)^{-1}$." What does this mean exactly?
And also what is a standard perturbation method for equations like  $y' - (y-1)^2 = \epsilon \frac{y^2}{x^2}$? 
I am only familiar with the case when $\epsilon$ only appears on the $y''$ term. Here I think we can not just plug in 
$y = y_0+ \epsilon y_1 + \epsilon^2 y_2+ \cdots $, since there are two possible boundary layers at $x = 0$ or $x=c$ when $y$ takes the non constant solution? 
And for what type of ODEs we can simply substitute in $y = y_0+ \epsilon y_1 + \epsilon^2 y_2+ \cdots $ to get the asymptotic?
 A: I cannot answer all of your questions, partially because I do not understand some of your questions (such as the "What does this mean exactly?" question, as this differential equation here is different from the one in the link).  I do  not understand why you say "Here I think we cannot just plug in [...]," since this procedure works perfectly fine in this problem.  However, in this answer, I shall show you that you can still use the standard perturbation theory to deal with this differential equation.
I shall write $y=f^\epsilon$ for the solution to
$$y'(x)-\big(y(x)-1\big)^2=\epsilon\,\frac{\big(y(x)\big)^2}{x^2}\text{ and }y(1)=1\,.\tag{*}$$
Note that $f^0$ is the constant function $f^0\equiv 1$.  Thus, for small $\epsilon$, we expect that $f^\epsilon$ is obtained from $f^0$ by via the perturbation theory.  That is, for some functions $f_0,f_1,f_2,\ldots$, we have
$$f^\epsilon=f_0+\epsilon\,f_1+\epsilon^2\,f_2+\ldots\,,$$
where $f_0=f^0$.  Observe that $f_k(1)=0$ for all $k=1,2,3,\ldots$.
The equation for $f_1$ is obtained by ignoring terms of oder $\epsilon^2$ or higher.   From (*) and from $f_0=f^0\equiv1$, we get
$$\epsilon\,f_1'(x)=\frac{\epsilon}{x^2}\text{ or }f_1'(x)=\frac{1}{x^2}\,.$$
Consequently, $$f_1(x)=\int_1^x\,\frac{1}{t^2}\,\text{d}t=1-\frac{1}{x}\,.$$
The equation for $f_2$ is obtained by ignoring terms of order $\epsilon^3$ or higher.  Since $f_0(x)=1$ and $f_1(x)=1-\frac{1}{x}$, we find that, from (*),
$$\left(\frac{\epsilon}{x^2}+\epsilon^2\,f'_2(x)\right)-\epsilon^2\left(1-\frac{1}{x}\right)^2=\frac{\epsilon}{x^2}\,\Biggl(1+2\,\epsilon\left(1-\frac1x\right)\Biggr)\,.$$
Therefore,
$$f_2'(x)-\left(1-\frac{1}{x}\right)^2=\frac{2}{x^2}\left(1-\frac1x\right)\,,\text{ or }f_2'(x)=1-\frac{2}{x}+\frac{3}{x^2}-\frac{2}{x^3}\,,$$
yielding
$$f_2(x)=\int_1^x\,\left(1-\frac{1}{t^2}\right)\,\text{d}t=(x-1)-2\,\ln(x)+3\,\left(1-\frac{1}{x}\right)-\left(1-\frac{1}{x^2}\right)\,.$$
You can continue so on and so forth, but the deeper you go, the more complicated the differential equation will be.  However, for small $\epsilon$ and for $x$ near $1$, you can see that the approximation
$$y(x)\approx 1+\epsilon\,\left(\frac{(x-1)}{x}\right)+\epsilon^2\,\Biggl((x-1)-2\,\ln(x)+3\,\left(1-\frac{1}{x}\right)-\left(1-\frac{1}{x^2}\right)\Biggr)$$
gives an estimate of the solution to (*) very well (check this with a numerical solver or something alike).  For $x$ further away from $1$ and for larger $\epsilon$, you are going to have to find more terms.
