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:


*

*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.

*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.

*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.

*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$.
 A: 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.

A: *

*Since this equation is nonlinear, there may also be distinguished limits where Y needs to be re-scaled as well as the independent variable x (so called corner layer). Also the boundary layer may be internal, so need xNew = (x-a)/eps^power for inner variable See J.D. Coles thin book (Peturbation methods, published in late sixties I think). Coles book has lots of inspired examples that date back to when this methodology was being developed for fluids problems.

*This particular equation is autonomous (no explicit x dependence), so very best way to solve it is by phase-planing. Can see how all the solutions change as A and B change with y(0) = A, y(1) = B

