diff --git a/doc/source/whatsnew/v1.0.0.rst b/doc/source/whatsnew/v1.0.0.rst
index 3f8d9d3916797..2a65c5c8eb458 100644
--- a/doc/source/whatsnew/v1.0.0.rst
+++ b/doc/source/whatsnew/v1.0.0.rst
@@ -550,6 +550,7 @@ Numeric
 - Bug in :class:`UInt64Index` precision loss while constructing from a list with values in the ``np.uint64`` range (:issue:`29526`)
 - Bug in :class:`NumericIndex` construction that caused indexing to fail when integers in the ``np.uint64`` range were used (:issue:`28023`)
 - Bug in :class:`NumericIndex` construction that caused :class:`UInt64Index` to be casted to :class:`Float64Index` when integers in the ``np.uint64`` range were used to index a :class:`DataFrame` (:issue:`28279`)
+- Bug in :meth:`Series.interpolate` when using method=`index` with an unsorted index, would previously return incorrect results. (:issue:`21037`)
 
 Conversion
 ^^^^^^^^^^
diff --git a/pandas/core/missing.py b/pandas/core/missing.py
index 744cde95cb668..9f4f7445509c9 100644
--- a/pandas/core/missing.py
+++ b/pandas/core/missing.py
@@ -277,7 +277,11 @@ def interpolate_1d(
                 inds = lib.maybe_convert_objects(inds)
         else:
             inds = xvalues
-        result[invalid] = np.interp(inds[invalid], inds[valid], yvalues[valid])
+        # np.interp requires sorted X values, #21037
+        indexer = np.argsort(inds[valid])
+        result[invalid] = np.interp(
+            inds[invalid], inds[valid][indexer], yvalues[valid][indexer]
+        )
         result[preserve_nans] = np.nan
         return result
 
diff --git a/pandas/tests/series/test_missing.py b/pandas/tests/series/test_missing.py
index 81bf1edbe86df..b1c9b67eb951e 100644
--- a/pandas/tests/series/test_missing.py
+++ b/pandas/tests/series/test_missing.py
@@ -1655,3 +1655,14 @@ def test_interpolate_timedelta_index(self, interp_methods_ind):
             pytest.skip(
                 "This interpolation method is not supported for Timedelta Index yet."
             )
+
+    @pytest.mark.parametrize(
+        "ascending, expected_values",
+        [(True, [1, 2, 3, 9, 10]), (False, [10, 9, 3, 2, 1])],
+    )
+    def test_interpolate_unsorted_index(self, ascending, expected_values):
+        # GH 21037
+        ts = pd.Series(data=[10, 9, np.nan, 2, 1], index=[10, 9, 3, 2, 1])
+        result = ts.sort_index(ascending=ascending).interpolate(method="index")
+        expected = pd.Series(data=expected_values, index=expected_values, dtype=float)
+        tm.assert_series_equal(result, expected)