Consider the following differential equation: $$ \frac{d^{2}y}{dx^{2}}=\frac{1}{2}\begin{cases} 1-e^{-\frac{y}{\epsilon}},\space\space\space x<0\\ e^{-\frac{1-y}{\epsilon}}-1,\space\space\space x>0 \end{cases} $$ subject to $\frac{dy}{dx}\to0$ as $y\to0,1$. I am interested in the behaviour of the solution in the limit $\epsilon\to0$. Solving this exactly leads to the implicit solution: $$ -\frac{x}{\epsilon}=\begin{cases} f\left(\frac{1}{2\epsilon},\frac{y}{\epsilon}\right), x<0\\ f\left(\frac{1-y}{\epsilon},\frac{1}{2\epsilon}\right), x>0 \end{cases} $$ where $$ f(y_{1},y_{0})=\int_{y_{0}}^{y_{1}}\frac{dy}{\sqrt{y+e^{-y}-1}} $$ I have been able to expand the solution to get $$ y\sim\begin{cases} A\epsilon e^{\frac{x+\sqrt{2}}{\sqrt{2\epsilon}}}, x<-\sqrt{2}\\ \frac{1}{4}(x+\sqrt{2})^{2}-\frac{x}{\sqrt{2}}\epsilon, -\sqrt{2}<x<0\\ 1-\frac{1}{4}(x-\sqrt{2})^{2}-\frac{x}{\sqrt{2}}\epsilon, 0<x<\sqrt{2}\\ 1-A\epsilon e^{-\frac{x-\sqrt{2}}{\sqrt{2\epsilon}}}, x>\sqrt{2} \end{cases} $$ where $A=e^{\frac{\int_{0}^{1}\left(\frac{1}{\sqrt{y+e^{-y}-1}}-\frac{\sqrt{2}}{y}\right)dy+\int_{1}^{\infty}\left(\frac{1}{\sqrt{y+e^{-y}-1}}-\frac{1}{\sqrt{y}}\right)dy-2}{\sqrt{2}}}\approx0.66237$.
My question is: how do I obtain this without knowledge of the exact solution? (the equation to which I want to find an asymptotic solution does not admit any form of exact solution). I can obtain the leading order solution by appropriately approximating the exponential terms in the differential equation (e.g. $e^{-\frac{y}{\epsilon}}\sim1$ where $y<O(\epsilon)$ and $e^{-\frac{y}{\epsilon}}\sim0$ where $y>O(\epsilon)$) to get: $$ y\sim\begin{cases} 0, x<-\sqrt{2}\\ \frac{1}{4}(x+\sqrt{2})^{2}, -\sqrt{2}<x<0\\ 1-\frac{1}{4}(x-\sqrt{2})^{2}, 0<x<\sqrt{2}\\ 0, x>\sqrt{2} \end{cases} $$ However, I'm struggling to find the first order correction. I'm trying a substitution of the form: $$ y(x)=y_{0}(x)+k\epsilon^{\alpha}e^{\frac{S(x)}{\beta}} $$ where $ y_{0}(x)$ is the above, piecewise quadratic, leading order solution and $S(x)$ is some $O(1)$ function to be determined, along with $k$, $\alpha$ and $\beta$, assuming $\alpha,\beta\geq0$ (but not both $0$). When I do this, I'm not able to solve for all of the unknowns and I don't see how the constant $A$ is going to be determined from the method. Any help is greatly appreciated.