Abstract Given a finite number of samples of a continuous set-valued function F, mapping an interval to compact subsets of the real line, we develop good approximations of F, which can be computed efficiently. In the first stage, we develop an efficient algorithm for computing an interpolant to $F$, inspired by the ‘metric polynomial interpolation’, which is based on the theory in Dyn et al. (2014, Approximation of Set-Valued Functions: Adaptation of Classical Approximation Operators. Imperial College Press). By this theory, a ‘metric polynomial interpolant’ is a collection of polynomial interpolants to all the ‘metric chains’ of the given samples of $F$. For set-valued functions whose graphs have nonempty interior, the collection of these ‘metric chains’ can be infinite. Our algorithm computes a small finite subset of ‘significant metric chains’, which is sufficient for approximating $F$. For the class of Lipschitz continuous functions with samples at the roots of the Chebyshev polynomials of the first kind, we prove that the error incurred by our computed interpolant decays with increasing number of interpolation points in the same rate as in the case of interpolation by the metric polynomial interpolant. This is also demonstrated by our numerical examples. For the class of set-valued functions whose graphs have smooth boundaries, we extend our algorithm to achieve a high-precision detection of the points of topology change, followed by a high-order approximation of the boundaries of the graph of F. We further discuss the case of set-valued functions whose graphs have ‘holes’ with Hölder-type singularities at the points of change of topology. To treat this case we apply some special approximation ideas near the singular points of the holes. We analyze the approximation order of the algorithm, including the error in approximating the points of change of topology, and show by several numerical examples the capability of obtaining high-order approximation of the holes.