diff --git a/lib/matplotlib/tests/test_wilkinson_locator.py b/lib/matplotlib/tests/test_wilkinson_locator.py
new file mode 100644
index 000000000000..eeb2f1b271dc
--- /dev/null
+++ b/lib/matplotlib/tests/test_wilkinson_locator.py
@@ -0,0 +1,188 @@
+import numpy as np
+import pytest
+from matplotlib.ticker import MaxNLocator
+from matplotlib.ticker import WilkinsonLocator
+
+
+def test_wilkinson_basic():
+    """Ticks should cover the full data range."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(3, 97)
+
+    assert len(ticks) >= 2
+    assert ticks[0] <= 3
+    assert ticks[-1] >= 97
+
+
+def test_wilkinson_vs_maxn():
+    """WilkinsonLocator should not produce more ticks than MaxNLocator."""
+    w = WilkinsonLocator(nbins=5).tick_values(3, 97)
+    m = MaxNLocator(nbins=5).tick_values(3, 97)
+
+    assert len(w) <= len(m)
+
+
+def test_wilkinson_reversed_input():
+    """Should handle vmin > vmax gracefully."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(97, 3)
+
+    assert len(ticks) >= 2
+    assert ticks[0] <= 3
+    assert ticks[-1] >= 97
+
+
+def test_wilkinson_equal_input():
+    """Should handle vmin == vmax without crashing."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(50, 50)
+
+    assert len(ticks) >= 1
+
+
+def test_wilkinson_negative_range():
+    """Should work correctly for negative ranges."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(-100, -10)
+
+    assert ticks[0] <= -100
+    assert ticks[-1] >= -10
+
+
+def test_wilkinson_cross_zero():
+    """Should work correctly when range crosses zero."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(-50, 50)
+
+    assert ticks[0] <= -50
+    assert ticks[-1] >= 50
+
+
+def test_wilkinson_small_range():
+    """Should work for very small ranges."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(0.001, 0.009)
+
+    assert len(ticks) >= 2
+    assert ticks[0] <= 0.001
+    assert ticks[-1] >= 0.009
+
+
+def test_wilkinson_large_range():
+    """Should work for very large ranges."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(0, 1_000_000)
+
+    assert len(ticks) >= 2
+    assert ticks[0] <= 0
+    assert ticks[-1] >= 1_000_000
+
+
+# ---------------- steps Parameter ---------------- #
+
+def test_wilkinson_custom_steps():
+    """Custom steps should be respected."""
+    loc = WilkinsonLocator(nbins=5, steps=[1, 2, 5, 10])
+    ticks = loc.tick_values(0, 100)
+
+    assert len(ticks) >= 2
+    assert ticks[0] <= 0
+    assert ticks[-1] >= 100
+
+
+def test_wilkinson_custom_steps_stored():
+    """Custom steps should be stored on the instance."""
+    custom = [1, 2, 5, 10]
+    loc = WilkinsonLocator(nbins=5, steps=custom)
+
+    assert loc.steps == custom
+
+
+def test_wilkinson_default_steps():
+    """Default steps should be [1, 2, 2.5, 5, 10]."""
+    loc = WilkinsonLocator(nbins=5)
+
+    assert loc.steps == [1, 2, 2.5, 5, 10]
+
+
+def test_wilkinson_steps_empty_raises():
+    """Empty steps list should raise ValueError."""
+    with pytest.raises(ValueError):
+        WilkinsonLocator(nbins=5, steps=[])
+
+
+def test_wilkinson_single_step():
+    """Single step value should still produce ticks."""
+    loc = WilkinsonLocator(nbins=5, steps=[1])
+    ticks = loc.tick_values(0, 100)
+
+    assert len(ticks) >= 2
+
+
+# ---------------- Scoring Fairness ---------------- #
+
+def test_wilkinson_coverage_dominates():
+    """
+    Coverage should dominate simplicity.
+    Both [0,25,50,75,100] and [0,20,40,60,80,100] are valid for (0,100).
+    The choice should depend on coverage/density, not just q niceness.
+    The key assertion is that the result covers the range well.
+    """
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(0, 100)
+
+    assert ticks[0] <= 0
+    assert ticks[-1] >= 100
+    assert len(ticks) >= 4
+
+
+def test_wilkinson_does_not_always_prefer_25_steps():
+    """
+    For a range like (0, 80), [0,20,40,60,80] (step=20) should score
+    higher than [0,25,50,75] (step=25) because it covers the range better.
+    """
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(0, 80)
+
+    # [0, 20, 40, 60, 80] covers exactly; [0, 25, 50, 75] misses 80
+    assert ticks[-1] >= 80
+
+
+def test_wilkinson_simplicity_not_sole_decider():
+    """
+    q=1 is always most 'simple', but shouldn't always win.
+    For (0, 100) with nbins=5, a step of 25 or 20 is better than step=1.
+    """
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(0, 100)
+
+    # If simplicity alone decided, we'd get 100 ticks of step=1
+    assert len(ticks) <= 10
+
+
+# ---------------- Tick Quality ---------------- #
+
+def test_wilkinson_ticks_are_sorted():
+    """Ticks should always be in ascending order."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(3, 97)
+
+    assert np.all(np.diff(ticks) > 0)
+
+
+def test_wilkinson_ticks_evenly_spaced():
+    """Ticks should be evenly spaced (uniform step)."""
+    loc = WilkinsonLocator(nbins=5)
+    ticks = loc.tick_values(0, 100)
+
+    diffs = np.diff(ticks)
+    assert np.allclose(diffs, diffs[0], rtol=1e-5)
+
+
+def test_wilkinson_nbins_respected():
+    """Number of ticks should stay close to nbins."""
+    for nbins in [3, 5, 8, 10]:
+        loc = WilkinsonLocator(nbins=nbins)
+        ticks = loc.tick_values(0, 100)
+        # Allow some flexibility but shouldn't be wildly off
+        assert len(ticks) <= nbins * 2 + 1
diff --git a/lib/matplotlib/ticker.py b/lib/matplotlib/ticker.py
index 83e13841677a..9a2096cb6d4d 100644
--- a/lib/matplotlib/ticker.py
+++ b/lib/matplotlib/ticker.py
@@ -154,7 +154,8 @@
            'Locator', 'IndexLocator', 'FixedLocator', 'NullLocator',
            'LinearLocator', 'LogLocator', 'AutoLocator',
            'MultipleLocator', 'MaxNLocator', 'AutoMinorLocator',
-           'SymmetricalLogLocator', 'AsinhLocator', 'LogitLocator')
+           'SymmetricalLogLocator', 'AsinhLocator', 'LogitLocator',
+           'WilkinsonLocator')
 
 
 class _DummyAxis:
@@ -3045,3 +3046,176 @@ def __call__(self):
     def tick_values(self, vmin, vmax):
         raise NotImplementedError(
             f"Cannot get tick locations for a {type(self).__name__}")
+
+
+class WilkinsonLocator(Locator):
+    """
+    Tick locator based on an Extended Wilkinson-style algorithm.
+    Balances simplicity, coverage, and density.
+
+    Parameters
+    ----------
+    nbins : int, optional
+        Target number of tick intervals. Default is 10.
+    steps : list of float, optional
+        Sequence of "nice" step values (mantissas) to consider.
+        Each value should be in the range [1, 10].
+        Default is [1, 2, 2.5, 5, 10], matching Wilkinson's original Q set.
+        Example: steps=[1, 2, 5, 10] for a simpler set.
+    """
+
+    def __init__(self, nbins=10, steps=None):
+        self._nbins = max(1, int(nbins))
+        if steps is None:
+            self.steps = [1, 2, 2.5, 5, 10]
+        else:
+            if not steps:
+                raise ValueError("steps must be a non-empty list")
+            self.steps = sorted([float(s) for s in steps])
+
+    def __call__(self):
+        vmin, vmax = self.axis.get_view_interval()
+        return self.tick_values(vmin, vmax)
+
+    def tick_values(self, vmin, vmax):
+    # Handle reversed inputs
+        if vmin > vmax:
+            vmin, vmax = vmax, vmin
+
+    # Handle singular case (vmin == vmax)
+        if vmin == vmax:
+            vmin -= 1
+            vmax += 1
+
+        locs = self._wilkinson(vmin, vmax)
+        return self.raise_if_exceeds(locs)
+
+    # ---------------- CORE ALGORITHM ---------------- #
+
+    def _wilkinson(self, vmin, vmax):
+        best_score = -np.inf
+        best_ticks = None
+
+        span = vmax - vmin
+        if span == 0:
+            return np.array([vmin])
+
+        k_min = int(np.floor(np.log10(span))) - 1
+        k_max = int(np.ceil(np.log10(span))) + 1
+
+        for q in self.steps:
+            for k in range(k_min, k_max + 1):
+                if abs(k) > 10:
+                    continue
+
+                step = q * (10 ** k)
+                if step <= 0:
+                    continue
+
+                lmin = np.floor(vmin / step) * step
+                lmax = np.ceil(vmax / step) * step
+
+                ticks = np.arange(lmin, lmax + 0.5 * step, step)
+                n = len(ticks)
+
+                # Reject too many or too few ticks
+                if n < 2 or n > self._nbins * 2:
+                    continue
+
+                score = self._score(vmin, vmax, ticks, q)
+
+                if score > best_score:
+                    best_score = score
+                    best_ticks = ticks
+
+        # Fallback
+        if best_ticks is None:
+            step = span / self._nbins
+            best_ticks = np.arange(vmin, vmax + step, step)
+
+
+        # Clamp to data range with small tolerance
+        best_ticks = best_ticks[
+            (best_ticks >= vmin - 1e-9) &
+            (best_ticks <= vmax + 1e-9)
+        ]
+
+        # Rebuild tick array using consistent step to avoid float drift
+        if len(best_ticks) >= 2:
+            step = best_ticks[1] - best_ticks[0]
+            start = np.floor(vmin / step) * step
+            end = np.ceil(vmax / step) * step
+            best_ticks = np.arange(start, end + 0.5 * step, step)
+
+        # Ensure at least 2 ticks
+        if len(best_ticks) < 2:
+            return np.array([vmin, vmax])
+
+        return best_ticks
+
+    # ---------------- SCORING ---------------- #
+
+    def _score(self, vmin, vmax, ticks, q):
+        """
+        Score a candidate tick sequence.
+
+        Weights are chosen so that coverage (how well ticks span the data)
+        and density (how close n is to nbins) dominate over simplicity
+        (niceness of q). This avoids always preferring q=2.5 (25-steps)
+        over q=2 (20-steps) regardless of context — the better choice
+        depends on coverage and density for the specific data range.
+        """
+        n = len(ticks)
+
+        simplicity = self._simplicity(q)
+        coverage = self._coverage(vmin, vmax, ticks)
+        density = self._density(n)
+
+        overflow = (
+            max(0, vmin - ticks[0]) +
+            max(0, ticks[-1] - vmax)
+        )
+
+        # Penalize exceeding target bin count
+        tick_penalty = max(0, n - self._nbins)
+
+        return (
+            0.15 * simplicity +   # Low weight: niceness of q alone shouldn't decide
+            0.50 * coverage +     # High weight: ticks must cover the data range well
+            0.25 * density -      # Medium weight: stay close to nbins
+            0.50 * tick_penalty -  # Penalize too many ticks
+            2.00 * overflow / (vmax - vmin)  # Penalize ticks outside data range
+        )
+
+    def _simplicity(self, q):
+        """
+        Score how 'nice' the step value q is.
+        Lower index in self.steps = simpler = higher score.
+        Returns a value in [0, 1].
+        """
+        n = len(self.steps)
+        if n == 1:
+            return 1.0
+        return 1 - self.steps.index(q) / (n - 1)
+
+    def _coverage(self, vmin, vmax, ticks):
+        """
+        Score how well the ticks cover [vmin, vmax].
+        Penalizes both under- and over-shooting the data range.
+        Returns a value close to 1 when ticks align tightly with data.
+        """
+        span = vmax - vmin
+        if span == 0:
+            return 1.0
+
+        return 1 - (
+            ((vmin - ticks[0]) ** 2 + (ticks[-1] - vmax) ** 2)
+            / (span ** 2)
+        )
+
+    def _density(self, n):
+        """
+        Score how close the number of ticks n is to the target nbins.
+        Returns 1 when n == nbins, decreasing as n moves away.
+        """
+        return max(0, 1 - abs(n - self._nbins) / self._nbins)