Approximate embeddings of the hyperbolic plane in $\Bbb R^3$ It is known that there is no isometric embedding of the hyperbolic plane into $\Bbb R^3$ that is $C^2$ or higher. (The Nash-Kuiper theorem guarantees the existence of an exact $C^1$ embedding.)
How does this situation change if we look at approximate embeddings, in which the embedded geodesic distance is always within some arbitrarily small $\epsilon$ of the true hyperbolic distance? Are $C^\infty$ approximate embeddings guaranteed to exist for any arbitrarily small $\epsilon$? If not, what is the smallest such $\epsilon$, and how does this change if we look at $C^2$ instead?
Formally, let's assume there is some map, call it $f$, which maps from $\Bbb H^2$ to $\Bbb R^3$, and that this mapping is suitably smooth (at least $C^2$, preferably $C^\infty$). Let's say that $d(a,b)$ is the distance between any two points in the hyperbolic plane, and $d_f(a,b)$ is the geodesic distance on the embedded surface in $\Bbb R^3$. Then an "$\epsilon$-approximate embedding" is one in which $|d_f(a,b) - d(a,b)| < \epsilon$, for all $a, b$.
Then the question is, what is the largest $n$ such that a $C^n$ approximate embedding of $\Bbb H^2 \to \Bbb R^3$ exists for all arbitrarily small $\epsilon$? And if such embeddings don't exist, are there any lower bounds on how small $\epsilon$ can be?
 A: There is no $\epsilon$-approximate embedding $f : \mathbb H^2 \to \mathbb R^3$, no matter what smoothness you pick nor what $\epsilon > 0$ you pick. 
In fact, forget smoothness altogether, there is simply no function $f : \mathbb H^2 \to \mathbb R^3$ such that $|d_f(a,b) - d(a,b)| < \epsilon$.
To see why, pick a large number $B$, and imagine, in $\mathbb H^2$, how many pairwise disjoint balls of radius $2\epsilon$ with can be packed into a single $B$ ball. The number is exponential in $B$, i.e. the maximum number of such balls has a lower bound of the form $c a^{B}$ where $c>0$ and $a>1$. Let's pick a set of balls that realizes this maximum, and let's enumerate their centers as $x_1,...,x_K$, $K > c a^{B}$. We have $d(x_k,x_{k'}) > 4 \epsilon$ for all $k \ne k' \in \{1,...,K\}$.
Now let's look at the radius $\epsilon$ balls $B(x_k,\epsilon) \subset \mathbb H^2$ and their images $f(B(x_k,\epsilon)) \subset \mathbb R^3$. Note that if $k \ne k'$ then 
$$d(f(x_k), f(x_{k'})) = d_f(x_k,x_{k'}) > d(x_k,x_{k'}) - \epsilon \ge 4\epsilon - \epsilon > 2\epsilon
$$
and therefore the $\epsilon$ balls around $f(x_k),f(x_{k'})$ are pairwise disjoint in $\mathbb R^3$. And the number of such balls has an exponential lower bound in $B$.
On the other hand, 
$$d(f(x_k),f(x_{k'})) \le d(x_k,x_{k'}) + \epsilon \le B + \epsilon
$$
The volume of the ball of radius $B+\epsilon$ is cubic in $B$. So we are fitting an exponential number of $\epsilon$-balls into a cubic volume, which for large $B$ is a contradiction.

Just as a comment, I believe that the Nash-Kuiper theorem is just about "infinitesmal" isometries, i.e. maps $f : \mathbb H^2 \to \mathbb R^3$ for which $D_x f : T_x \mathbb H^2 \to T_{f(x)} \mathbb R^3$ is an isometric embedding of inner product spaces, for each $x$. Their theorem is ruling out any such map which is $C^2$, regardless of whether or not it satisfies a distance lower bound of the form you have hypothesized, i.e. $d_f(x,y) > d(x,y)-\epsilon$. My "easy" proof uses that lower bound strongly.
A: The noncompactness of $\mathbb{H}^2$ combined with the rather stringent notion of "$\epsilon$-approximate embedding" you consider make your question rather challenging. That is why I will not properly answer your question here, but I still wish to bring up some clues, some 'food for thoughts', that for any $\epsilon > 0$, there might exist a $C^{\infty}$-embedding $f : \mathbb{H}^2 \to \mathbb{R}^3$ such that the distance induced by the Riemannian metric $f^*g_{euclidean}$ is uniformly close to the hyperbolic metric.
Below, given a Riemannian manifolds $(M, g)$ and $(N,h)$, a $C^1$-map $f : (M, g) \to (N,g)$ will be said to be isometric if $f^*h = g$.
Context. Given a closed (i.e. compact boundarily) surface $\Sigma$ and a $C^2$-immersion $f : \Sigma \to \mathbb{R}^3$ (where 3-space is equipped with the Euclidean Riemannian metric), the curvature of $f^*g_{euc}$ is well-defined and in fact equals the determinant of the (well-defined!) shape operator (that is Gauss' theorema egrigium). By compactness, this determinant has to be strictly positive at least at one point. Consequently, there is no isometric $C^2$-immersion of a strictly negatively curved closed surface $(\Sigma, g)$ into $\mathbb{R}^3$. This fact was generalised by Hilbert in 1901: given a complete Riemannian surface $(\Sigma, g)$ of constant strictly negative curvature, there is no isometric sufficiently regular immersion $f : (\Sigma, g) \to (\mathbb{R}^3, g_{euc})$. This fact was itself later generalised by Efimov: the same conclusion holds whenever the curvature of $g$ is smaller some strictly negative constant. (See T.K. Milnor's lecture for details.) Interestingly, given any Riemannian surface $(\Sigma, g)$, Nash-Kuiper theorem states that there exists a $C^1$-isometric embedding $f : (\Sigma, g) \to (U, g_{euc})$ where $U \subset \mathbb{R}^3$ is any open set (the shape operator of $f$ can't necessarily be defined in this case, so the theorema egrigium is no obstruction).
Nash-Kuiper's approximation method. An isometric map $f : (\Sigma, g) \to (\mathbb{R}^3, g_{euc})$ is a solution to, a 'root' of the equation $f^*g_{euc} - g = 0$. The vague idea behind Nash and Kuiper result is the following. They begin with an approximate $C^{\infty}$-solution $f_0 : (\Sigma, g) \to (\mathbb{R}^3, g_{euc})$ to $f^*g_{euc}-g=0$ and then apply a generalisation of Newton's method to get a $C^1$-convergent sequence of $C^{\infty}$-maps $f_j$ which are ever better approximations to this equation. In the process, they lose all control on second and higher-orders derivatives. This approximation method was generalized by Gromov to yield his convex-integration technique, which is briefly explained at the end of Eliashberg-Mishachev Introduction to the $h$-principle (this part of the book can be read independently from the rest of the book; a proof of the Nash-Kuiper theorem is given for compact manifolds).
Of course, one needs to specify in what sense these functions are 'approximations'. In Eliashberg-Mishachev, they consider quasi-isometries, namely maps $f$ such that $C^{-1} g \le f^*g_{euc} \le C g$ for some constant $C \ge 1$. In the compact case they consider, this is equivalent to your '$\epsilon$-approximation' notion, but it is more general than your notion when the source manifold is non-compact. Perhaps the proof can be adapted in order to replace $C$ by a rapidly decreasing function. In any case, this suggests that what you seek for is true, since it is for closed surfaces.
More direct approaches. These $h$-principle techniques are quite general, but that makes them difficult to visualize even in specific situations as that of $(\Sigma, g) = (\mathbb{H}^2, g_{hyp})$. One could perhaps argue more directly.
(1) As another small clue of the possibility of "$\epsilon$-approximations" of the hyperbolic plane into 3-space, consider the upper-half plane model of $\mathbb{H}^2$ and more specifically the subsets $H = \{(x,y) \in \mathbb{R}^2 | y > 1 \}$ and $H' = \{(x,y) \in H | x \in [0, 2\pi) \}$. $H$ is a "large" portion of $\mathbb{H}^2$ to the extend that it is the interior of a horocycle. The map $p : H \to H' : (x,y) \mapsto (x \mbox{ mod } 2\pi, y)$ is an isometry. Now, there exist an embedding $f' : H' \to \mathbb{R}^3$ whose image is the pseudosphere. The pseudosphere admits a normal unit vector field $\vec{N}$. Given an appropriate diffeomorphism $g : \mathbb{R} \to (-\delta, \delta)$, the map $f : H \to \mathbb{R}^3 : (x,y) \mapsto f'(p(x,y)) + g(x)\vec{N}(f'(p(x,y)))$ is an $\epsilon$-approximate isometry.
(2) On pages 49-50 of his book Three-dimensional geometry and topology, Thurston describes (in an exercise) how to physically construct approximate (in a quasi-isometric sense) portions of the hyperbolic plane out of paper. This construction gained the interest of Henderson and Taimina, who wrote a lot about it. Taimina found out how to crochet this construction out of yarn. Taimina produced in this way several very interesting and colorful figures; see for instance this webpage.
The most mathematical and detailed description I managed to find about this construction is the following text by Henderson and Taimina. It seems to be implied there and in other texts from them that if one had an infinite amount of time, one could produce in this way arbitrarily good quasi-isometric embedding of the whole of $\mathbb{H}^2$.
The representation of the hyperbolic space they use is not quite the upper half-space model, but rather $(\mathbb{R}^2, g)$ defined by pulling back $g_{hyp}$ by the map $h : \mathbb{R}^2 \to H : (x,y) \mapsto (x, e^y)$. (I assume here, that $g_{hyp}$ has curvature -1.) Note that (or read the first answer here) that the horocycles $y=y_0$ have geodesic curvature 1, like the circle of radius 1 inside the Euclidean 2-plane. This suggests that the "fattened horocycles" $\mathbb{R} \times [y_0, y_0 + \delta]$ can be approximated by (the universal covers of) the annuli $r \in [1, 1 + \delta'(y_0, \delta)]$ in the Euclidean 2-plane. One then glues together those annuli to get an approximate representation of the hyperbolic plane. (To me, this suggests that in Thurston's construction, if one wants to have an $\epsilon$-approximate isometry and not merely a quasi-isometric set, the width of the annuli should vary. However, in Henderson-Taimina's text I mentioned above, they say "You can experiment with different ratios BUT not in the same model. You will get a hyperbolic plane ONLY if you will be increasing the number of stitches in the same ratio all the time." So I don't know what to think.)
