Given a set A of n points in Rd with weight function w : A → R>0, the Fermat distance function is φ(x):= Pa∈A w(a)∥x−a∥. A classic problem in facility location dating back to 1643, is to find the Fermat point x∗, the point that minimizes the function φ. We consider the problem of computing a point xe∗ that is an ε-approximation of x∗ in the sense that ∥xe∗ − x∗∥ < ε. The algorithmic literature has so far used a different notion based on ε-approximation of the value φ(x∗). We devise a certified subdivision algorithm for computing xe∗, enhanced by Newton operator techniques. We also revisit the classic Weiszfeld-Kuhn iteration scheme for x∗, turning it into an ε-approximate Fermat point algorithm. Our second problem is the certified construction of ε-isotopic approximations of n-ellipses. These are the level sets φ−1(r) for r > φ(x∗) and d = 2. Finally, all our planar (d = 2) algorithms are implemented in order to experimentally evaluate them, using both synthetic as well as real world datasets. These experiments show the practicality of our techniques.