TL;DR: Your options are limited 1) go brute force adaptive for accurate and expensive solution 2) use numerical diffusion for a less accurate but stable solution or (my favorite) 3) leverage the fact that this is a singular perturbation problem and solve two inexpensive inner/outer problems and let matched asymptotics do its magic!
If you really must obtain a uniform numerical solution to the problem, there is really not much you could do beyond adaptive mesh refinement. You are facing a singular perturbation problem that develops a boundary layer of thickness $\delta = \mathcal{O}(\sqrt{\epsilon})$ near the boundary. This boundary layer separates the inner and outer solutions.
To be more precise, consider the left boundary, where $x = \mathcal{O}(\delta)$. In this region, it is possible to rescale the equation, by using stretched coordinate $\eta = x/\delta$, for the inner solution:
$$
-\Delta u_i + u_i = 1
$$
subjected to the boundary conditions $u(0) = 0$ and $u_i(\eta\rightarrow\infty) = u_o(x\rightarrow0)$. For the outer solution, $u_o$, one assumes that $x = \mathcal{O}(1)$ and thus the dominant balance is between $u$ and $1$ which simply gives $u_0 = 1$. With the outer solution at hand, you can now solve the regularized inner solution with ease -- in this case even analytically.
This is in fact the technique that was (and still is) very popular for solving laminar boundary layer problems in fluid mechanics back in the day. In fact if you look at the Navier-Stokes equations, at high Reynolds numbers, you are effectively facing a singular perturbation problem, which like the one you mentioned here, develops a boundary layer (fun fact: the terms "boundary layer" in perturbation analysis actually comes from the fluid boundary layer problem I just described).
So, back in the day where computational resources were very limited, the typical thing to do was to separate the outer and inner problems. For the outer, one would first solve the potential flow (analogous to $u_0 = 1$ here) and for the inner, well the boundary layer equations! All this without the needing fancy adaptive methods if you were to tackle the whole problem at once and in a single computational domain.