diff --git a/docs/changes.rst b/docs/changes.rst
index 95335a5..60d2e6d 100644
--- a/docs/changes.rst
+++ b/docs/changes.rst
@@ -5,12 +5,12 @@ PyDL Changelog
 1.0.0 (unreleased)
 ------------------
 
-* Update testing infrastructure (PR `#67`_).
-* Update links to documentation (PR `#64`_).
-* Remove support for 2-dimensional B-spline, which is better supported in
+* Deprecate support for 2-dimensional B-spline, which is better supported in
   other packages such as PypeIt_ or specreduce_ (PR `#73`_).
 * Use specutils_ as an informal template for Python packaging boilerplate (PR `#73`_).
 * Change Harris Geospatial to NV5 Geospatial.
+* Update testing infrastructure (PR `#67`_).
+* Update links to documentation (PR `#64`_).
 
 .. _`#64`: https://github.com/weaverba137/pydl/pull/64
 .. _`#67`: https://github.com/weaverba137/pydl/pull/67
diff --git a/pydl/pydlutils/bspline.py b/pydl/pydlutils/bspline.py
index 8263698..9544d86 100644
--- a/pydl/pydlutils/bspline.py
+++ b/pydl/pydlutils/bspline.py
@@ -542,6 +542,9 @@ def iterfit(xdata, ydata, invvar=None, upper=5, lower=5, x2=None,
             maxiter=10, **kwargs):
     """Iteratively fit a B-spline set to data, with rejection.
 
+    Additional keyword parameters are passed to
+    :class:`~pydl.pydlutils.bspline.bspline`.
+
     Parameters
     ----------
     xdata : :class:`numpy.ndarray`
@@ -626,6 +629,8 @@ def iterfit(xdata, ydata, invvar=None, upper=5, lower=5, x2=None,
     ywork = ydata[xsort]
     invwork = invvar[xsort]
     if x2 is not None:
+        warn('2D bspline fits may be buggy and will be fully removed in the future.',
+             DeprecationWarning)
         x2work = x2[xsort]
     else:
         x2work = None
diff --git a/pydl/pydlutils/tests/test_bspline.py b/pydl/pydlutils/tests/test_bspline.py
index d70d0a2..9da44dd 100644
--- a/pydl/pydlutils/tests/test_bspline.py
+++ b/pydl/pydlutils/tests/test_bspline.py
@@ -42,6 +42,38 @@ def test_iterfit():
     # print(yfit)
     # pylab.plot(x, y, 'k-', x, yfit, 'r-')
 
+def test_iterfit_2d():
+    y0 = np.array([0.661984, 0.134913, 0.0410350, 0.940134, 0.411034,
+                   0.484675, 0.169943, 0.325046, 0.269194, 0.552381,
+                   0.797177, 0.971658, 0.251765, 0.531675, 0.854556,
+                   0.411237, 0.694380, 0.499562, 0.437242, 0.362451,
+                   0.343206, 0.524099, 0.158634, 0.728597, 0.198340,
+                   0.571210, 0.477527, 0.962797, 0.973921, 0.413651,
+                   0.736380, 0.516366, 0.104283, 0.675993, 0.467053,
+                   0.230112, 0.866994, 0.469885, 0.964392, 0.541084,
+                   0.332984, 0.581252, 0.422322, 0.872555, 0.803636,
+                   0.520998, 0.918942, 0.241564, 0.169263, 0.686649,
+                   0.708284, 0.707858, 0.00113957, 0.827920, 0.845985,
+                   0.416961, 0.553842, 0.526549, 0.501051, 0.337514,
+                   0.700873, 0.152816, 0.762935, 0.650039, 0.483321,
+                   0.708600, 0.410033, 0.507671, 0.596956, 0.177692,
+                   0.498112, 0.422037, 0.788333, 0.856578, 0.941245,
+                   0.432411, 0.356469, 0.341916, 0.0331059, 0.641100,
+                   0.690452, 0.168667, 0.915178, 0.158406, 0.701508,
+                   0.841774, 0.434161, 0.153123, 0.420066, 0.0499331,
+                   0.947241, 0.0768818, 0.410540, 0.843788, 0.0640255,
+                   0.513463, 0.511104, 0.680434, 0.762480, 0.0563867])
+    assert y0.size == 100
+    y = smooth(y0, 10)
+    assert y.size == 100
+    x = np.arange(y0.size, dtype='d')
+    x2 = x.copy()
+    with pytest.warns(DeprecationWarning) as w:
+        sset, outmask = iterfit(x, y, x2=x2, maxiter=0, nord=3, bkspace=10)
+    assert len(w) > 0
+    assert str(w[0].message) == '2D bspline fits may be buggy and will be fully removed in the future.'
+    assert sset.npoly == 1
+    assert sset.funcname == 'legendre'
 
 def test_cholesky_band():
     ab = np.array([[8., 9., 10., 11., 12., 13., 14., 0., 0., 0.],
diff --git a/pydl/tests/test_pydl.py b/pydl/tests/test_pydl.py
index bf8c7a9..60a4850 100644
--- a/pydl/tests/test_pydl.py
+++ b/pydl/tests/test_pydl.py
@@ -262,7 +262,11 @@ def test_smooth():
 def test_uniq():
     items = np.array([1, 2, 3, 1, 5, 6, 1, 7, 3, 2, 5, 9, 11, 1])
     items_sorted = np.sort(items)
-    items_argsorted = np.argsort(items)
+    #
+    # On some systems, kind='quicksort', which is the default, gives different
+    # results depending on numpy version.
+    #
+    items_argsorted = np.argsort(items, kind='stable')
     #
     # Test pre-sorted array.
     #