The outer solution $y_{\text{outer}}(x) = -1$ is fine, but it turns out that in this problem the width of the boundary layer isn't $\approx \epsilon$ as you have assumed with your substitution $x = \epsilon X$. To determine the correct width, instead substitute $x = \epsilon^{\alpha} X$, where $\alpha$ is a constant to be determined. The equation
$$
\epsilon y'' + xy' + \epsilon y = 0
$$
then becomes
$$
\epsilon^{1-2\alpha} \frac{d^2 Y}{dX^2} + X \frac{dY}{dX} + \epsilon Y = 0.
$$
The third term, $\epsilon Y$, being order $\epsilon$, is certainly smaller than the second term, $X \frac{dY}{dX}$, which is order $1$. Therefore the third term must not be involved in the dominant balance. This only leaves the first two terms, which are balanced when $\alpha = 1/2$. The resulting equation is then, to leading order,
$$
\frac{d^2 Y}{dX^2} + X \frac{dY}{dX} = 0.
$$
After multiplying by the integrating factor $e^{X^2/2}$ we can write the equation as
$$
\frac{d}{dX} \left(e^{X^2/2} \frac{dY}{dX}\right) = 0,
$$
and integrating yields
$$
\frac{dY}{dX} = A e^{-X^2/2}.
$$
Integrating once more we see that
$$
\newcommand{erf}{\operatorname{erf}}
Y(X) = A\erf\left(\frac{X}{\sqrt{2}}\right) + B,
$$
where $\erf$ is the error function. To satisfy the boundary condition $y(0) = 1$ we must then take $B = 1$. To match this with the outer solution we need
$$
\lim_{X \to \infty} Y(X) = y_{\text{outer}}(0) = -1,
$$
which is satisfied when $A = -2$.
We now have the inner and outer solutions, namely
$$
y_{\text{outer}}(x) = -1, \\
y_{\text{inner}}(x) = 1 - 2\erf\left(\frac{x}{\sqrt{2\epsilon}}\right).
$$
The matched solution is therefore
$$
y(x) \sim y_{\text{inner}}(x) + y_{\text{outer}}(x) - y_{\text{outer}}(0) = 1 - 2\erf\left(\frac{x}{\sqrt{2\epsilon}}\right).
$$
Here's a plot which shows the actual solution in blue and this approximate solution in purple for $\epsilon = 1/50$.
Equivalent forms of the equation
With the substitution $y=Au(Bt)$ scaling both the domain and the range, the equation reduces to
$$
u''+\frac{\pi^2+ε}{B^2}u+\frac{A}{B^2}u^2=0.
$$
The scales can now be chosen in a way that produces the unperturbed equation with $ε=0$ using $A=B^2=1+\frac{ε}{\pi^2}$, or a form closer to the standard harmonic oscillator
$$
x''+x+δx^2=0
$$
with $B^2=\pi^2+ε$, $A=δB^2$.
Citing a previously performed perturbation analysis
This last form, after switching the sign in $x$, was explored in a previous topic. Fixing the amplitude $R$, the solution close to the circle of radius $R$ is given by the parametrization
$$
x\sqrt{1+\frac23δx}=-R\cos\phi, ~~ \dot x=R\sin\phi,\\
\implies x=-R\cos\phi-δ\frac{R^2}3\cos^2\phi+O(δ^2),\\
\left[\small\begin{aligned}
\implies 1&=(1+\frac{2δR}3\cos\phi+...)\dot\phi,\\
t&=\phi+\frac{2δR}3\sin\phi+...\iff \phi=t-\frac{2δR}3\sin t+...\\
x(t)&=-R\left[\cos t+\sin t\left(\frac{2δR}3\sin t\right)+...\right]-δ\frac{R^2}3\cos^2t+O(δ^2)\\
&=-R\cos t-δ\frac{R^2}3[\cos^2t+2\sin^2 t]+O(δ^2)
\end{aligned}\right]
$$
for the half-period of interest $\phi\in[0,\pi]$. Taking the derivative of the first equation one can derive a first-order ODE for the angle. Doing some series expansion one gets an expression for the half-period $\phi(T)=\pi$ as
$$
T=\pi\left(1+\frac{5}{12}(δR)^2+\frac{385}{576}(δR)^4+O(δ^6)\right).
$$
Per the parametrization we want $B=T$, which in the leading order gives
$$
(\piδR)^2=\frac65ε.
$$
This suggests $δ=\sqrtε$, $R=\frac{\sqrt{6/5}}{\pi}$.
Fall back to perturbation series for numerics
Now set $x=-R\cos t+δx_1$ to result in the equation
$$
x_1''+x_1+x^2=0,
$$
which in first order has the initial condition $x_1(0)=-\frac{R^2}3$ and the solution
$$
x_1=\frac{R^2}6(\cos(2t)-3).
$$
This can be used in an approximation $x=x_0+δx_1+δ^2v$ or just as initial guess in the boundary value problem for $v$ in $x=x_0+δv$.
Numerical experiments
# global constants and invariant functions
R2 = 1.2/np.pi**2
R = R2**0.5
def x0(t): return -R*np.cos(t)
def dx0(t): return R*np.sin(t)
def x1(t): return R2/6*(np.cos(2*t)-3)
def dx1(t): return -R2/3*np.sin(2*t)
# Variables depending on the perturbation parameter
eps = 1e-8
dta = eps**0.5
B2 = np.pi**2+eps
B, A = B2**0.5, dta*B2
dvb = -dx0(B)/dta
# Boundary Value Problem
def ode_v(t,v): xt= x0(t)+dta*v[0]; return [v[1], -v[0]-xt**2]
def bc_v(v0,vb): return [v0[1],vb[1]-dvb]
# initial guess
t = np.linspace(0,B,5)
v = [ x1(t), dx1(t) ]
# solve and print (error) message
res = solve_bvp(ode_v,bc_v,t,v, tol=1e-4)
print(res.message,res.y[0,0], res.y[0,-1],-R2/3)
This converges well and reliably. Plotting after scaling back for the original problem gives in the phase plane
Best Answer
As you approach $y=1$, or better $y=1-ε$, the next term to become "active" is the third one. So try $y=\frac1{1+εu}$ with $u>1+O(ε)$ to get $$ -εu'=ε^2u+1-\frac{1+εu}{u}=\frac{ε^2u^2+(1-ε)u-1}{u}\tag1 $$ Thus approximately $$ -εu'=1-\frac1u\iff u'+\frac{u'}{u-1}=-\frac1ε \\ u(t)+\ln(u(t)-1)=-\frac{t}ε+d \\ (u(t)-1)e^{u(t)-1}=De^{-\frac{t}ε}\iff u(t)=1+W(De^{-\frac{t}ε}) ,~~~ y(t)=\frac1{1+ε+εW(De^{-\frac{t}ε})}\tag2 $$ using the Lambert-W function
This is a rather short segment of width $ε$ around the $t_0$ given in the question.
Then even closer to the asymptote, set $u=1+δv$ and insert in (1) $$ -εδv'=ε^2(1+δv)+1-\frac{1}{1+δv}-ε=ε^2(1+δv)+\frac{δv}{1+δv}-ε. $$ This balances on the right for $δ=ε$ to give in the leading order terms $$ -εv'=v-1\implies v(t)=1+Fe^{-tε},~~~ y(t)=\frac1{1+ε+ε^2+ε^2Fe^{-tε}} $$
Now select the constants so that all fits together...