The usual 2-sphere exists naturally in $\mathbb R^3$, and in general the usual definition of $S^n$ is as a particular subset of $\mathbb R^{n+1}$ with the induced metric. In that case, the identity map is a locally metric-preserving embedding into $\mathbb R^2$, but it doesn't preserve the global distance. To wit, two diametrically opposed points have distance $2$ in $\mathbb R^3$ but distance $\pi$ along geodesics in the sphere itself.
Thus, the natural embedding works as an isometry when we view the two spaces as Riemannian manifolds, but not when we consider them directly as metric spaces. It appears that both kinds of maps can be called "isometric embeddings", but nonetheless they are different concepts.
In differential geometry, a hyperbolic plane is a complete, simply-connected Riemannian $2$-manifold equipped with a metric of constant negative (Gaussian) curvature. By scaling the metric, we may assume the curvature is $-1$; often when people say hyperbolic plane they're assuming $K \equiv -1$. The rest of this answer assumes $K \equiv -1$.
It turns out that any two hyperbolic planes of curvature $-1$ are isometric as Riemannian manifolds. In my experience, when people speak of the hyperbolic plane, they refer to the equivalence class of hyperbolic planes in the space of Riemannian manifolds up to isometry.
It's possible that a hyperbolic model in the sense of the question is what I've called a hyperbolic plane above. There are a number of standard models. The first three below naturally sit inside Minkowski three-space, the real Cartesian three-space equipped with the metric $dx^{2} + dy^{2} - dz^{2}$.
The hyperboloid model, the upper sheet of a two-sheeted hyperboloid, $x^{2} + y^{2} - z^{2} = -1$, $z > 0$. Geodesics ("hyperbolic lines") in this model are intersections of the hyperboloid with planes through the origin. (This surface has non-constant _positive curvature with respect to the Euclidean ambient metric, but acquires constant negative curvature from the ambient Minkowski metric.)
The conformal disk model (or Poincaré model), the open unit disk in the plane $z = 0$ with the metric $\frac{4(dx^{2} + dy^{2})}{(1 - (x^{2} + y^{2}))^{2}}$ induced by projection from the point $(0, 0, -1)$, see diagram. Because the metric is a scalar multiple of the Euclidean metric $dx^{2} + dy^{2}$, hyperbolic angles coincide with Euclidean angles. Geodesics turn out to be arcs of circles meeting the unit circle at right angles. At first glance this metric appears flat; see below.
The affine disk model (or Klein-Beltrami model), the open unit disk in the plane $z = 1$ with the metric induced by projection from the origin. Geodesics in this model are Euclidean segments. If the open unit disk and projection point are translated one unit up the $z$-axis in the preceding model, we obtain the affine model.
The conformal upper half plane model, induced from the conformal disk by a fractional linear transformation. The metric is $\frac{dx^{2} + dy^{2}}{y^{2}}$, and geodesics are vertical lines and semicircles centered on the $x$-axis.
Curvature is a measure of angular defect per unit area in geodesic triangles. Qualitatively, a surface has negative curvature if the total interior angle of a geodesic triangle (enclosing a topological disk, which is automatic in the hyperbolic plane) is smaller than $\pi$. Given that hyperbolic and Euclidean angles coincide in the conformal models, and that geodesics are arcs of circles, it's visually plausible that the conformal disk and half-plane models have negative curvature.
If we represent a region of the hyperbolic plane as a surface in Euclidean three-space, we obtain a surface having "saddle-like" geometry at each point. The best-known examples are surfaces of rotation, especially the pseudosphere, and surfaces with helicoidal symmetry such as Dini's surface.
A hyperbolic circle of hyperbolic radius has hyperbolic circumference $2\pi \sinh r = \pi(e^{r} - e^{-r})$. Particularly, circumference grows exponentially with radius. The theorem of Hilbert mentioned asserts that the entire hyperbolic plane cannot be represented, even allowing self-intersection. Qualitatively (thinking of crocheted models that start like a saddle and grow outward radially), a hyperbolic disk is "floppy", and the larger the radius, the more area must be accommodated in a Euclidean ball whose Euclidean radius grows like the hyperbolic radius of the disk; "Euclidean space just can't keep up".
As noted in the comments, however, a hyperbolic disk of arbitrarily large hyperbolic radius can be isometrically embedded in Euclidean three-space. One method is to use Dini's surface, taking the edge of the disk to lie on the edge of the surface, and taking the center to lie "exponentially far away along the horn", so that the disk is tightly wrapped around the "axis".
Coda: Hyperbolic geometry also overlaps synthetic geometry. Euclid's parallel postulate is replaced by another, equivalent to "Given a line $\ell$ and a point $p$ not on $\ell$, there exist two lines through $p$ that do not meet $\ell$." Pat Ryan's Euclidean and Non-Euclidean Geometry is an accessible treatment, including the hyperboloid model and the necessary geometry of Minkowski space, and requiring little more than one semester of linear algebra.
Best Answer
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.)