diff --git a/numpy/polynomial/polyutils.py b/numpy/polynomial/polyutils.py
index 5c65e03c2f4b..25d50837a3c6 100644
--- a/numpy/polynomial/polyutils.py
+++ b/numpy/polynomial/polyutils.py
@@ -289,8 +289,8 @@ def mapparms(old, new) :
     Parameters
     ----------
     old, new : array_like
-        Each domain must (successfully) convert to a 1-d array containing
-        precisely two values.
+        Domains. Each domain must (successfully) convert to a 1-d array
+        containing precisely two values.
 
     Returns
     -------
@@ -330,13 +330,14 @@ def mapdomain(x, old, new) :
     """
     Apply linear map to input points.
 
-    The linear map ``offset + scale*x`` that maps `old` to `new` is applied
-    to the points `x`.
+    The linear map ``offset + scale*x`` that maps the domain `old` to
+    the domain `new` is applied to the points `x`.
 
     Parameters
     ----------
     x : array_like
-        Points to be mapped.
+        Points to be mapped. If `x` is a subtype of ndarray the subtype
+        will be preserved.
     old, new : array_like
         The two domains that determine the map.  Each must (successfully)
         convert to 1-d arrays containing precisely two values.
@@ -388,6 +389,6 @@ def mapdomain(x, old, new) :
     array([-1.0+1.j , -0.6+0.6j, -0.2+0.2j,  0.2-0.2j,  0.6-0.6j,  1.0-1.j ])
 
     """
-    [x] = as_series([x], trim=False)
+    x = np.asanyarray(x)
     off, scl = mapparms(old, new)
     return off + scl*x
diff --git a/numpy/polynomial/tests/test_polyutils.py b/numpy/polynomial/tests/test_polyutils.py
index 86f2a5b9b3b0..14bf8bb78984 100644
--- a/numpy/polynomial/tests/test_polyutils.py
+++ b/numpy/polynomial/tests/test_polyutils.py
@@ -67,9 +67,25 @@ def test_mapdomain(self) :
         dom1 = [0 - 1j, 2 + 1j]
         dom2 = [-2, 2]
         tgt = dom2
-        res = pu.mapdomain(dom1, dom1, dom2)
+        x = dom1
+        res = pu.mapdomain(x, dom1, dom2)
         assert_almost_equal(res, tgt)
 
+        # test for multidimensional arrays
+        dom1 = [0,4]
+        dom2 = [1,3]
+        tgt = np.array([dom2, dom2])
+        x = np.array([dom1, dom1])
+        res = pu.mapdomain(x, dom1, dom2)
+        assert_almost_equal(res, tgt)
+
+        # test that subtypes are preserved.
+        dom1 = [0,4]
+        dom2 = [1,3]
+        x = np.matrix([dom1, dom1])
+        res = pu.mapdomain(x, dom1, dom2)
+        assert_(isinstance(res, np.matrix))
+
     def test_mapparms(self) :
         # test for real values
         dom1 = [0,4]