I have found the following singular perturbation problem,
$\epsilon u_{xx} + |u_x|u_x + u = 0$, $x>0$;
with initial conditions, $u(0) = \epsilon^2$, $u_x(0) = 0$, where $0 < \epsilon \ll 1$. The question is to find an approximate solution $\forall x > 0$.
I am rather new with these methods and I have never found a problem like this (I have just worked with more or less basic problems where you can easily apply multiple scales method or WKB, for example), so I do not have much idea about how to proceed. Any hint?
Written as a two-dimensional dynamical system, the ODE yields \begin{align} u' &= v,\\ \epsilon v' &= -u - v |v|. \tag{1} \end{align} From a dynamical systems point of view, this system is a bit unusual since the vector field is $C^1$ but not $C^2$. However, that shouldn't be a problem per se. Since you're interested in an initial condition very close to the origin, in particular $(u(0),v(0)) = (\epsilon^2,0)$, it is useful to look at the stability of the origin. The eigenvalues of the linearisation of $(1)$ around $(0,0)$ are $\pm \frac{i}{\sqrt{\epsilon}}$, which means the origin is a center. Therefore, we can't invote the Grobman-Hartman theorem directly, but we see that the nonlinear term $-v|v|$ is attractive. Therefore, exactly as @David suggested, we expect an oscillatory solution with slowly decreasing amplitude. Moreover, the oscillation frequency is moderately fast, as it scales as $\frac{1}{\sqrt{\epsilon}}$ (see the eigenvalues of the origin).
From the initial condition $u(0)=\epsilon^2$, we see that it is useful to introduce a rescaled variable $u = \epsilon^2 y$. If we wouldn't do that, and look for solutions of $\mathcal{O}(1)$ or $\mathcal{O}(\epsilon)$, the initial condition would just be $u(0)= 0$ to leading order, yielding the trivial solution $u = 0$.
Anyway, in terms of $y$, the ODE is \begin{equation} \epsilon y_{xx} + \epsilon^2 y_x |y_x| + y = 0,\quad y(0)=1,\;y_x(0)=0. \tag{2} \end{equation} Introducing the rescaled variable $\xi = \frac{x}{\sqrt{\epsilon}}$, we get \begin{equation} y_{\xi \xi} + \epsilon y_\xi |y_\xi| + y = 0,\quad y(0)=1,\;y_\xi(0)=0. \tag{3} \end{equation} Now this is an equation we can work with! Clearly, the leading order solution to $(3)$ is $\cos \xi$. But how to get an expression for the slowly decreasing amplitude?
The best way to approach this, I think, is to use the method of averaging. For more information, see F. Verhulst, Nonlinear Differential Equations and Dynamical Systems (2nd ed.), Springer, 2006, chapter 11. To be brief, our ODE $(3)$ is of the form \begin{equation} y_{\xi\xi} + y = \epsilon f(y,y_\xi) = - \epsilon y_\xi |y_\xi|. \end{equation} The solutions to the unperturbed ODE are $\alpha \cos \xi$ and $\beta \sin \xi$. Suppose we write the solution to the full ODE as \begin{equation} y(\xi) = \alpha(\xi) \cos \xi + \beta(\xi) \sin \xi, \end{equation} such that \begin{equation} y'(\xi) = -\alpha(\xi) \sin\xi + \beta(\xi) \cos \xi \end{equation} (we can always make this choice by choosing suitable $\alpha$ and $\beta$). Then, we obtain \begin{align} -\alpha_\xi \sin\xi + \beta_\xi \cos\xi &= - \epsilon \left((\beta(\xi) \cos \xi -\alpha(\xi) \sin\xi)|\beta(\xi) \cos \xi -\alpha(\xi) \sin\xi|\right),\\ \alpha_\xi \cos \xi + \beta_\xi \sin\xi &= 0, \end{align} yielding \begin{align} \alpha_\xi &= \epsilon \sin \xi \left((\beta(\xi) \cos \xi -\alpha(\xi) \sin\xi)|\beta(\xi) \cos \xi -\alpha(\xi) \sin\xi|\right),\\ \beta_\xi &= -\epsilon \cos \xi \left((\beta(\xi) \cos \xi -\alpha(\xi) \sin\xi)|\beta(\xi) \cos \xi -\alpha(\xi) \sin\xi|\right). \end{align} Now, we average the right hand side over one period in $\xi$, leaving $\alpha$ and $\beta$ fixed. This yields a system for the averaged values of $\alpha$ and $\beta$, denoted as $\alpha_a$ and $\beta_a$, \begin{align} \alpha_{a,\xi} &= - \epsilon \frac{8}{3} \alpha_a \sqrt{\alpha_a^2+\beta_a^2},\\ \beta_{a,\xi} &= - \epsilon \frac{8}{3} \beta_a \sqrt{\alpha_a^2+\beta_a^2}. \end{align} We can solve this system easily using polar coordinates, i.e. $\alpha_a = r_a \cos \theta_a$ and $\beta_a = r_a \sin \theta_a$, yielding \begin{align} r_{a,\xi} &= -\epsilon \frac{8}{3} r_a^3,\\ \theta_{a,\xi} &= 0. \end{align} Thus, the general leading order solution to $(3)$ is \begin{equation} y(\xi) = \frac{\cos (\xi-\xi_0)}{\sqrt{\frac{16}{3}\epsilon\xi + r_0}}. \end{equation} The initial conditions fix $\xi_0 = 0$ and $r_0 = 1$. In terms of the original spatial variable $x$ and the original function $u$, we thus obtain \begin{equation} u(x) = \epsilon^2\frac{\cos \frac{x}{\sqrt{\epsilon}}}{\sqrt{\frac{16}{3}\sqrt{\epsilon}x+1}} \end{equation} as a leading order, multiple scale approximation of the solution to the original ODE.