diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index 5c595bcf668d..90011de99ffc 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -1777,27 +1777,165 @@ def chebroots(c):
     return r
 
 
-def chebinterpolate(func, deg, args=()):
+def __dct(x, n=None):
+    """
+    Discrete Cosine Transform type II
+                      N-1
+           y[k] = 2 * sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
+                      n=0
+    Examples
+    --------
+    >>> x = np.arange(5)
+    >>> x1 = np.arange(6)
+    >>> y = x + 1j*x
+    >>> y1 = x1 + 1j*x1
+
+    >>> np.allclose(__dct(x), [20. , -9.95959314,  0. , -0.89805595,  0. ])
+    True
+
+    >>> np.allclose(__dct(x1),
+    ...             [30. , -14.41953704, 0. , -1.41421356, 0. , -0.27740142])
+    True
+
+    >>> np.allclose(__dct(y), [20.   +20.j, -9.95959314 -9.95959314j,
+    ...                         0.    +0.j, -0.89805595 -0.89805595j,
+    ...                         0.    +0.j        ])
+    True
+
+    >>> np.allclose(__dct(y1), [ 30. +30.j , -14.41953704-14.41953704j,
+    ...                           0.  +0.j ,  -1.41421356 -1.41421356j,
+    ...                           0.  +0.j ,  -0.27740142 -0.27740142j])
+    True
+
+    >>> np.allclose(x, __idct(__dct(x)))
+    True
+
+    >>> np.allclose(x, __dct(__idct(x)))
+    True
+
+    References
+    ----------
+    .. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
+           Makhoul, `IEEE Transactions on acoustics, speech and signal
+           processing` vol. 28(1), pp. 27-34,
+           :doi:`10.1109/TASSP.1980.1163351` (1980).
+    .. [2] Wikipedia, "Discrete cosine transform",
+           https://en.wikipedia.org/wiki/Discrete_cosine_transform
+    """
+    fft = np.fft.fft
+    x = np.atleast_1d(x)
+
+    if n is None:
+        n = x.shape[-1]
+
+    if x.shape[-1] < n:
+        n_shape = x.shape[:-1] + (n - x.shape[-1],)
+        xx = np.hstack((x, np.zeros(n_shape)))
+    else:
+        xx = x[..., :n]
+
+    real_x = np.all(np.isreal(xx))
+    if (real_x and (np.remainder(n, 2) == 0)):
+        xp = 2 * fft(np.hstack((xx[..., ::2], xx[..., ::-2])))
+    else:
+        xp = fft(np.hstack((xx, xx[..., ::-1])))[..., :n]
+
+    w = np.exp(-1j * np.arange(n) * np.pi / (2 * n))
+
+    y = xp * w
+
+    if real_x:
+        return y.real
+    return y
+
+
+def __idct(x, n=None):
+    """
+    Inverse Discrete Cosine Transform type II
+                N-1
+    y[k] = 1/N sum w[n]*x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
+               n=0
+    w[0] = 1/2
+    w[n] = 1 for n>0
+
+    Examples
+    --------
+    >>> x = np.arange(5)
+    >>> np.allclose(x, __idct(__dct(x)))
+    True
+
+    >>> np.allclose(x, __dct(__idct(x)))
+    True
+
+    References
+    ----------
+    .. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
+           Makhoul, `IEEE Transactions on acoustics, speech and signal
+           processing` vol. 28(1), pp. 27-34,
+           :doi:`10.1109/TASSP.1980.1163351` (1980).
+    .. [2] Wikipedia, "Discrete cosine transform",
+           https://en.wikipedia.org/wiki/Discrete_cosine_transform
+    """
+
+    ifft = np.fft.ifft
+    x = np.atleast_1d(x)
+
+    if n is None:
+        n = x.shape[-1]
+
+    w = np.exp(1j * np.arange(n) * np.pi / (2 * n))
+
+    if x.shape[-1] < n:
+        n_shape = x.shape[:-1] + (n - x.shape[-1],)
+        xx = np.hstack((x, np.zeros(n_shape))) * w
+    else:
+        xx = x[..., :n] * w
+
+    real_x = np.all(np.isreal(x))
+    if (real_x and (np.remainder(n, 2) == 0)):
+        xx[..., 0] = xx[..., 0] * 0.5
+        yp = ifft(xx)
+        y = np.zeros(xx.shape, dtype=complex)
+        y[..., ::2] = yp[..., :n // 2]
+        y[..., ::-2] = yp[..., n // 2::]
+    else:
+        yp = ifft(np.hstack((xx,
+                             np.zeros_like(xx[..., 0]),
+                             np.conj(xx[..., :0:-1]))))
+        y = yp[..., :n]
+
+    if real_x:
+        return y.real
+    return y
+
+
+def chebinterpolate(func, deg, args=(), domain=None):
     """Interpolate a function at the Chebyshev points of the first kind.
 
-    Returns the Chebyshev series that interpolates `func` at the Chebyshev
-    points of the first kind in the interval [-1, 1]. The interpolating
-    series tends to a minmax approximation to `func` with increasing `deg`
-    if the function is continuous in the interval.
+    Returns the Chebyshev series that interpolates the multi-dimensional
+    function `func` at the Chebyshev points of the first kind scaled and
+    shifted to the `domain`. The interpolating series tends to a minmax
+    approximation of `func` with increasing `deg` if the function is
+    continuous in the input domain.
 
     .. versionadded:: 1.14.0
 
     Parameters
     ----------
     func : function
-        The function to be approximated. It must be a function of a single
-        variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are
-        extra arguments passed in the `args` parameter.
-    deg : int
-        Degree of the interpolating polynomial
+        The function to be approximated. It must be a function of the form
+        ``f(x1, x2,..., xn, a, b, c...)``, where ``x1, x2,..., xn`` are the
+        n-dimensional coordinates and ``a, b, c...`` are extra arguments
+        passed in the `args` parameter.
+    deg : list of int, size ndim
+        Degree of the interpolating polynomial used for each dimension.
     args : tuple, optional
         Extra arguments to be used in the function call. Default is no extra
         arguments.
+    domain : array-like, shape ndim x 2, optional
+        defines the domain ``[a1, b1] x [a2, b2] x ... x [an, bn]`` over
+        which `func` is interpolated. The default is None, in which case
+        the domain is [-1,1] * len(deg))
 
     Returns
     -------
@@ -1808,10 +1946,17 @@ def chebinterpolate(func, deg, args=()):
     Examples
     --------
     >>> import numpy.polynomial.chebyshev as C
-    >>> C.chebfromfunction(lambda x: np.tanh(x) + 0.5, 8)
-    array([  5.00000000e-01,   8.11675684e-01,  -9.86864911e-17,
-            -5.42457905e-02,  -2.71387850e-16,   4.51658839e-03,
-             2.46716228e-17,  -3.79694221e-04,  -3.26899002e-16])
+    >>> c1 = C.chebinterpolate(lambda x: np.tanh(x) + 0.5, 8)
+    >>> x = C.chebpts1(9)
+    >>> np.allclose(C.chebval(x, c1), np.tanh(x)+0.5))
+    True
+
+    Two-dimensional interpolation to machine accuracy
+    >>> x1,x2 = np.meshgrid(x,x)
+    >>> c2 = C.chebinterpolate(lambda x,y: np.tanh(x+y) + 0.5, deg=(8, 8))
+    >>> np.allclose(C.chebval2d(x1, x2, c2), np.tanh(x1+x2)+0.5))
+    True
+
 
     Notes
     -----
@@ -1824,24 +1969,41 @@ def chebinterpolate(func, deg, args=()):
     zero. For instance, if the function is even then the coefficients of the
     terms of odd degree in the result can be set to zero.
 
+    References
+    ----------
+    .. [1] 'A Survey of Methods of Computing Minimax and Near-Minimax
+           Polynomial Approximations for Functions of a Single Independent
+           Variable' by W. Fraser, `Journal of the ACM` (JACM),
+           vol. 12(3), pp. 295-314, (1965)
+    .. [2] MathWorld, 'Chebyshev Approximation Formula',
+           https://mathworld.wolfram.com/ChebyshevApproximationFormula.html
+    .. [3] Wikipedia, 'Chebyshev nodes',
+           http://en.wikipedia.org/wiki/Chebyshev_nodes
+
     """
-    deg = np.asarray(deg)
+    deg = np.atleast_1d(deg)
 
     # check arguments.
-    if deg.ndim > 0 or deg.dtype.kind not in 'iu' or deg.size == 0:
+    if deg.ndim != 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
         raise TypeError("deg must be an int")
-    if deg < 0:
+    if np.any(deg < 0):
         raise ValueError("expected deg >= 0")
-
     order = deg + 1
-    xcheb = chebpts1(order)
-    yfunc = func(xcheb, *args)
-    m = chebvander(xcheb, deg)
-    c = np.dot(m.T, yfunc)
-    c[0] /= order
-    c[1:] /= 0.5*order
+    ndim = len(order)
 
-    return c
+    if domain is None:
+        domain = chebdomain.tolist() * ndim
+    domain = np.atleast_2d(domain).reshape((-1, 2))
+    x_i = [pu.mapdomain(chebpts1(o), chebdomain, d)
+           for o, d in zip(order, domain)]
+    c = func(*np.meshgrid(*x_i), *args)
+
+    for i in range(ndim):
+        c = __dct(c[..., ::-1])
+        c[..., 0] /= 2.
+        if i < ndim - 1:
+            c = np.moveaxis(c, source=-1, destination=0)
+    return c / np.product(order)
 
 
 def chebgauss(deg):
@@ -2060,13 +2222,12 @@ def interpolate(cls, func, deg, domain=None, args=()):
 
         Notes
         -----
-        See `numpy.polynomial.chebfromfunction` for more details.
+        See `numpy.polynomial.chebyshev.chebinterpolate` for more details.
 
         """
         if domain is None:
             domain = cls.domain
-        xfunc = lambda x: func(pu.mapdomain(x, cls.window, domain), *args)
-        coef = chebinterpolate(xfunc, deg)
+        coef = chebinterpolate(func, deg, args=args, domain=domain)
         return cls(coef, domain=domain)
 
     # Virtual properties
diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py
index 2f54bebfdb27..24c4906e6c65 100644
--- a/numpy/polynomial/tests/test_chebyshev.py
+++ b/numpy/polynomial/tests/test_chebyshev.py
@@ -506,6 +506,20 @@ def powx(x, p):
                 c = cheb.chebinterpolate(powx, deg, (p,))
                 assert_almost_equal(cheb.chebval(x, c), powx(x, p), decimal=12)
 
+    def test_2d_approximation(self):
+
+        def powxy(x, y, p):
+            return (x+y)**p
+
+        x = np.linspace(-1, 1, 10)
+        x1, x2 = np.meshgrid(x, x)
+        for deg in range(0, 10):
+            for p in range(0, deg + 1):
+                c = cheb.chebinterpolate(powxy, (deg, deg), (p,))
+                assert_almost_equal(cheb.chebval2d(x1, x2, c),
+                                    powxy(x1, x2, p),
+                                    decimal=12)
+
 
 class TestCompanion: