Find a composite solution to the following problem

I'm trying to solve the following exercise:

Find a composite solution to the following problem: $$ \epsilon y'' + y(y' + 3) = 0 \text{ for }0<x<1, \text{ where }y(0) = 1, \,y(1) = 1 $$ where $\epsilon <<1.$

What I've tried: I've just learned about matched asymptotic expansions and I'm pretty sure that my teacher wants me to use those to solve this exercise. The procedure works as follows:

1. Find the outer solution. This solution is often found by assuming that the solution can be expanded in powers of $\epsilon$. Because this solution often has only one arbitrary constant it will not be able to satisfy both the boundary value conditions.
2. Find the inner solution. Assume that there exists a boundary layer at one of the boundaries. Introduce a boundary-layer coordinate and use an expansion for the boundary-layer solution.
3. Matching. Since the inner (corresponding to the boundary-layer) and outer expansions are approximations of the same function we expect them to be the same in the region between the inner and outer layers.
4. Composite expansion. Use the inner and outer solution to find a solution that will work over the entire interval. This can be done by adding the two solutions and subtracting corresponding to the region where they are equal.

I found the outer solution by assuming that the solution $y$ can expanded in powers of $\epsilon$ as follows: $$ y \sim y_0(x) + \epsilon y_1(x) + \ldots \tag{1} $$ If we substitute $(1)$ into the problem equation we get $$ \epsilon(y''_0 + \epsilon y''_1 + \ldots) + (y_0 + \epsilon y_1 + \ldots)(y'_0 + y'_1 + \ldots + 3) = 0 $$ We can find $y_0$ by looking at the order one terms:

$\mathcal{O}(1):$ $$ y_0(y'_0 + 3) = 0 $$ so that $y_0$ or $y_0 = c_1 - 3x$. This is the first moment where I don't exactly know how to proceed: I would have expected one solution.

To find the inner solution or boundary-layer solution we assume that there is a boundary layer at $x = 0$ and introduce a boundary layer coordinate: $$ \bar{x} = \dfrac{x}{\epsilon^\alpha} $$ From this change of variables and the chain rule, we have that $$ \dfrac{d}{dx} = \dfrac{d\bar{x}}{dx}\dfrac{d}{d\bar{x}} = \dfrac{1}{\epsilon^\alpha}\dfrac{d}{d\bar{x}} $$ If we now let $Y(\bar{x})$ denote the solution when using this boundary-layer coordinate, the problem equation transforms to $$ \epsilon^{1 - 2\alpha}\dfrac{d^2}{d\bar{x}^2}(Y_0 + \epsilon Y_1 + \ldots) + (Y_0 + \epsilon Y_1 + \ldots)\dfrac{1}{\epsilon^\alpha}\dfrac{d}{d\bar{x}}(Y_0 + \epsilon Y_1 + \ldots) + 3(Y_0 + \epsilon Y_1 + \ldots) $$ In order to determine $Y_0$, I want to balance the terms in this equation. I tried balancing the first and second term so that the third term becomes higher order. To balance the first and second term we need $1 - 2\alpha = -\alpha$ so that $\alpha = 1$. In this case the first to terms are $\mathcal{O}(\frac{1}{\epsilon})$ and the third term is $\mathcal{O}(1)$ so that the requirement that the third term is higher order is satisfied. If $\alpha = 1$ we find:

$\mathcal{O}(\frac{1}{\epsilon})$:

\begin{align} Y''_0 + Y_0Y'_0 = 0\\ Y_0(0) = 1 \end{align} If I enter this equation in wolfram I get a pretty long solution and I don't think it's supposed to be like that. It feels as if I'm doing something wrong.

Question: I am doing this correctly? What should I be doing differently? I don't know how to proceed from here since I can't find a reasonable solution for $Y_0$.

The solution $y=0$ is incompatible with $y(0)=1$ or $y(1)=1$, thus you get as outer solution $y(x)=1-3x$ or $y(x)=4-3x$.

The balancing looks good. The resulting equation $Y''+YY'=0$ is independent of the location of the boundary layer. It can be directly integrated to $$ Y'+\frac12Y^2=C $$ This only has bounded solutions for $X\to\infty$ if $C=\frac12a^2$ is positive, $a \ge 1$. Then $$ \ln|Y+a|-\ln|Y-a|=aX+b $$

boundary layer at $x=0$

With $y_{outer}(x)=4-3x$ the compatibility condition is $4=y_{outer}(0)=Y(+\infty)=a$, and $1=Y(0)\implies b=\ln\frac53$. Thus the inner solution $y_{inner}=Y(x/ϵ)$ can be computed as $$ \frac{4+Y}{4-Y}=\frac53e^{4X} %3Y+12=20e-5eY, Y(3+5e)=20e-12 \implies Y(X)=4\,\frac{5e^{4X}-3}{5e^{4X}+3} $$

In composition, the first order approximation reads as $$ y(x)=y_{outer}(x)-y_{outer}(0)+y_{inner}(x)=4\,\frac{5-3e^{-4x/ϵ}}{5+3e^{-4x/ϵ}}-3x $$

Let's see how that compares to a numerical solution (using scipy.integrate.solve_bvp)

boundary layer at $x=1$

The outer solution is $y=1-3x$. The inner solution for $x=1+ϵX$ has now to satisfy $-2=y_{outer}(1)=Y(-\infty)=-a$, and $Y(0)=1\implies b=\ln 3$. This gives the inner solution $y_{inner}(x)=Y((x-1)/ϵ)$ as $Y(X)=2\,\frac{3e^{2X}-1}{3e^{2X}+1}$ so that in combination the inner-outer approximation is $$ y(x)=y_{outer}(x)-y_{outer}(1)+y_{inner}(x)=2\,\frac{3e^{-2(1-x)/ϵ}-1}{3e^{-2(1-x)/ϵ}+1}+3(1-x) $$

In the numerical experiments it is much harder to achieve convergence towards these solutions, the solver prefers the boundary layer at $x=0$. The first graph is from a not completely converged solution, the other solutions are for values of $ϵ$ that are a factor of $10$ smaller than in the first plot.