Find the composite solution to this problem

I want to find a composite solution to the boundary value problem: $$\epsilon y'' - y' + y^2 = 1, \text{ for }0 where $$\epsilon\ll 1$$.

My approach: I know that I can find a composite solution in four steps:

1. Find an outer solution.
2. Find a boundary-layer solution.
3. Apply matching so that the outer solution and the boundary-layer solution both approximate the same function correctly.
4. Find the composite solution, by adding the outer solution and the boundary-layer solution and subtracting the part where they are equal.

Step 1: Outer solution:

I assume that the solution has an expansion in powers of $$\epsilon$$. So that $$y\sim y_0 + \epsilon^\alpha y_1 + \epsilon^\beta y_2 + \ldots$$ with $$0 < \alpha < \beta <\ldots$$ If we substitute this into the equation we get $$\epsilon(y_0 + \epsilon^\alpha y_1 + \ldots )'' - (y_0 + \epsilon^\alpha y_1 + \ldots)' + (y_0 + \epsilon^\alpha y_1 + \ldots )^2 = 1$$ If we only look at the $$\mathcal{O}(1)$$ terms we get $$-y_0' + y_0^2 = 1$$ so that the outer solution becomes $$y_0 = \dfrac{1 - e^{2c_1 + 2x}}{e^{2c_1 + 2x} + 1}$$ Step 2: Boundary layer solution: Let's assume that there is a boundary layer at $$x = 0$$. I introduce the boundary layer coordinate $$\bar{x} = \dfrac{x}{\epsilon^\alpha} \Leftrightarrow \dfrac{d}{dx} = \dfrac{1}{\epsilon^\alpha}\dfrac{d}{d\bar{x}}, \dfrac{d^2}{dx^2} = \dfrac{1}{\epsilon^{2\alpha}}\dfrac{d^2}{d\bar{x}^2}$$ If we let $$Y(\bar{x})$$ denote the solution of the problem when using this boundary layer coordinate, the original equation becomes $$\epsilon^{1 - 2\alpha}\dfrac{d^2}{d\bar{x}^2}(Y_0 + \epsilon^\gamma Y_1 + \ldots) - \epsilon^{-\alpha}\dfrac{d}{d\bar{x}}(Y_0 + \epsilon^\gamma Y_1 + \ldots) + (Y_0 + \epsilon^\gamma Y_1 + \ldots ) = 1$$ To balance this equation we need to look at the different terms. We already used the second and third term in part 1 for the outer solution. We can try balancing the first and the second term so that the third term becomes higher order. We need $$1 - 2\alpha = -\alpha\Leftrightarrow \alpha = 1$$. With $$\alpha = 1$$ the equation with the boundary layer coordinate becomes 

$$\dfrac{1}{\epsilon}\dfrac{d^2}{d\bar{x}^2}(Y_0 + \epsilon^\gamma Y_1 + \ldots ) - \dfrac{1}{\epsilon}\dfrac{d}{d\bar{x}}(Y_0 + \epsilon^\gamma Y_1 + \ldots) + (Y_0 + \epsilon^\gamma Y_1 + \ldots)^2 = 1$$ If we now look at the order $$\mathcal{O}\big(\dfrac{1}{\epsilon}\big)$$ terms we get: $$Y_0''(\bar{x}) - Y_0'(\bar{x}) = 1, \, Y_0(0) = 1/3$$ so that $$Y_0(\bar{x}) = c_1e^{\bar{x}} + c_2 - \bar{x}$$. We need $$Y_0(0) = 1/3$$ so we have $$c_1 + c_2 = 1/3$$. If there is a boundary layer at $$x = 0$$ then we need the outer solution to satisfy the boundary condition at $$x = 1$$ so that we have $$\dfrac{1 - e^{2c_1 + 2}}{e^{2c_1 + 2} + 1} = 1$$. This last expression can be rewritten to \begin{align} 1 - e^{2c_1+2} = e^{2c_1 + 2} + 1\\ \Leftrightarrow -2e^{2c_1 + 2} = 0 \end{align} Since there is no finite $$c_1$$ for which this expression holds I should probably look if the outer solution should instead satisfy the boundary condition at $$x = 0$$. In that case we need to find $$c_1$$ such that $$-2e^{2c_1 + 2} = 1/3$$ but this results in a complex value for $$c_1$$ and I don't think that I should end up with such a solution.

Question: What am I doing wrong? How can I find a correct outer solution and a correct boundary layer solution so that I can start with the matching process?

The outer solution with $$y(1)=1$$ is $$y(x)=1$$, which you can find as one of the stationary points $$y_*=\pm 1$$ of $$y'=y^2-1$$. In general you should get $$y(x)=\frac{1-C_1e^{2x}}{1+C_1e^{2x}}.$$
The inner equation is the same independent of the basis point, as the equation is autonomous. You missed to multiply the right side with $$ϵ$$, which then gives the equation as $$Y''(X)-Y'(X)=0$$. The solution then is $$Y(X)=C_2e^X+C_3.$$ This is not bounded in direction $$X\to+\infty$$ so that the boundary layer can only be at the right boundary.
To satisfy the left boundary condition in the outer solution you need $$C_1=\frac12$$. Thus you need a jump from $$C_3=\frac{2-e^2}{2+e^2}$$ to $$C_2+C_3=1$$ on the right boundary.
The combined approximation is then $$y(x)=\frac{2-e^{2x}}{2+e^{2x}}+\left(1-\frac{2-e^2}{2+e^2}\right)e^{(x-1)/ϵ}$$ 