What is the realization of this surface in Euclidean space? I am wondering how to realize this abstract surface as a subset of Euclidean space.
The surface as a point set is the 2 dimenional Euclidean plane minus the origin.
The metric is given by declaring the following 2 vector fields to be an orthonormal frame:
$$
e_1 = x\ ∂_x - y\ ∂_y
\qquad\text{and}\qquad
e_2 = x^2\ ∂_x + e^{−xy}\ ∂_y
$$
I think that this surface has constant negative curvature equal to -1 but I do not think it is a Poincare disk minus a point because the origin is infinitely far away from any point on the x or y axes.
BTW: Here is the curvature calculation.
the Lie bracket, [e1,e2], equals e2 so the covariant derivative can be defined using the following formulas:
∇e1e1 = 0
∇e1e2 = 0
∇e2e1 = -e2
∇e2e2 = e1
Extend these formulas to arbitrary vector fields by linearity and the Leibniz rule.
The connection 1 form is ω(e2) = -1, ω(e1) = 0.
Since this connection is torsion free the exterior derivative of ω easily computes and is the volume element of the metric. Therefore the Gauss curvature is -1.
It may be helpful to note that the rectangular hyperpolas y = k/x are geodesics as are the rays from the origin along the axes.
 A: 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).
