In trying to understand the possible realizations (i.e., isometric embeddings) of your Riemannian surface into Euclidean $3$-space, the best thing to do is to first figure out how the domain can be immersed into a complete model of a simply-connected surface of constant negative curvature (i.e., the Poincaré half-plane), and then study the problem by considering the known (local) isometric embeddings of (pieces of) the Poincaré plane into Euclidean $3$-space. This two-step procedure is likely to yield more insight than trying to do it all in one step.
In fact, you'd be better off thinking of this in terms of Lie groups, the connected non-abelian Lie group of dimension $2$ in particular. This is because the geometry on the domain in question is completely determined by the pair of basis vector fields $e_1$ and $e_2$ that satisfy the equation $[e_1,e_2] = e_2$. They form a $2$-dimensional Lie algebra and thus define the action of a local Lie group on the domain. In your particular case, they aren't complete, but since the domain is simply connected, that's enough to guarantee that the domain you have can be immersed into the standard domain with this example.
To compute this 'explicitly', here is what you do: Let $f(t) = t + e^{-t}$ and consider the pair of differential equations
$$
\frac{da}{a} = \frac{(f(xy)-xy)\ dx - x^2\ dy}{x\ f(xy)}
\qquad\text{and}\qquad
\frac{db}{1} = a\frac{y\ dx + x\ dy}{x\ f(xy)}
$$
Because the $1$-form on the right hand side of the first equation is closed (no matter what smooth function $f$ of one variable you put there), as long as a connected domain $D$ where $x f(xy) \not=0$ is simply connected (which holds in your case), there is no problem seeing that $a(x,y)>0$ exists on $D$ and is determined up to multiplication by a positive constant. [In fact, one has $a(x,y) = x\ g(xy)$, where $g(t)$ satisfies the ode $f(t)g'(t) + g(t) = 0$.] Once $a$ is determined satisfying the first equation, the right hand side of the second equation is also seen to be closed, so, again, assuming that $D$ is connected and simply connected, there will exist a $b(x,y)$ satisfying that equation, and it will be unique up to an additive constant. [In fact, one has $b(x,y) = h(xy)$, where $h(t)$ satisfies the differential equation $h'(t) = g(t)/f(t)$.]
So suppose $a$ and $b$ satisfying the above equations have been found. Then the mapping $F:D\to \mathbb{R}^+\times \mathbb{R}$ defined by
$$
F(x,y) = \bigl(a(x,y),b(x,y)\bigr) = \bigl(\ x\ g(xy),\ h(xy)\ \bigr)
$$
will immerse $D$ into the right $ab$-half-plane isometrically with respect to your given metric on $D$ and the Poincaré metric $ds^2 = (da^2+db^2)/a^2$. It will also carry the vector field $e_1$ to $a\partial_a$ and the vector field $e_2$ to $a\partial_b$.
Your question amounts to asking what the map $F$ looks like when $f(t) = t + e^{-t}$ and $D$ is either the right half-plane $x>0$ or the left half-plane $x<0$. (This is because $f(t)>0$ for all $t$.) Obviously, the map $F$ carries the hyperbolae $xy=c$ to horizonal lines in the right $ab$-half-plane, but the remainder of the description depends on the ranges of the functions $g$ and $h$ that are derived from $f$. In the particular case at hand, when $f(t) = t + e^{-t}$, it is not hard to show that $g$ and $h$ are each bounded above and below. (You have to do a little estimating to show this because $g$ and $h$ are not elementary functions.) In particular, this implies that the image of $F$ is a horizontal strip in the right $ab$-half-plane, i.e., it is described by inequalities $a>0$ and $b_- < b < b_+$ for some constants $b_- < b_+$. All such strips are equivalent under isometries of the Poincaré metric. Geometrically, it represents the region between two geodesics that limit to a common point on the ideal boundary (which happens to be at infinity in this presentation).
