diff --git a/doc/changes/devel/13063.bugfix.rst b/doc/changes/devel/13063.bugfix.rst
new file mode 100644
index 00000000000..76eba2032a1
--- /dev/null
+++ b/doc/changes/devel/13063.bugfix.rst
@@ -0,0 +1 @@
+Fix bug in the colorbars created by :func:`mne.viz.plot_evoked_topomap` by `Santeri Ruuskanen`_.
\ No newline at end of file
diff --git a/examples/visualization/evoked_topomap.py b/examples/visualization/evoked_topomap.py
index 83d1916c6f9..53b7a60dbba 100644
--- a/examples/visualization/evoked_topomap.py
+++ b/examples/visualization/evoked_topomap.py
@@ -5,8 +5,8 @@
 Plotting topographic maps of evoked data
 ========================================
 
-Load evoked data and plot topomaps for selected time points using multiple
-additional options.
+Load evoked data and plot topomaps for selected time points using
+multiple additional options.
 """
 # Authors: Christian Brodbeck <christianbrodbeck@nyu.edu>
 #          Tal Linzen <linzen@nyu.edu>
diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index 10ec5459e02..b047de4ea32 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -27,6 +27,7 @@
     _clean_names,
     _is_numeric,
     _pl,
+    _time_mask,
     _to_rgb,
     _validate_type,
     fill_doc,
@@ -1988,10 +1989,18 @@ def plot_evoked_joint(
     contours = topomap_args.get("contours", 6)
     ch_type = ch_types.pop()  # set should only contain one element
     # Since the data has all the ch_types, we get the limits from the plot.
-    vmin, vmax = ts_ax.get_ylim()
+    vmin, vmax = (None, None)
     norm = ch_type == "grad"
     vmin = 0 if norm else vmin
-    vmin, vmax = _setup_vmin_vmax(evoked.data, vmin, vmax, norm)
+    time_idx = [
+        np.where(
+            _time_mask(evoked.times, tmin=t, tmax=None, sfreq=evoked.info["sfreq"])
+        )[0][0]
+        for t in times_sec
+    ]
+    scalings = topomap_args["scalings"] if "scalings" in topomap_args else None
+    scaling = _handle_default("scalings", scalings)[ch_type]
+    vmin, vmax = _setup_vmin_vmax(evoked.data[:, time_idx] * scaling, vmin, vmax, norm)
     if not isinstance(contours, list | np.ndarray):
         locator, contours = _set_contour_locator(vmin, vmax, contours)
     else:
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index d83698acbb1..d2c28c0a937 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -2116,6 +2116,22 @@ def plot_evoked_topomap(
     :ref:`gridspec <matplotlib:arranging_axes>` interface to adjust the colorbar
     size yourself.
 
+    The defaults for ``contours`` and ``vlim`` are handled as follows:
+
+    * When neither ``vlim`` nor a list of ``contours`` is passed, MNE sets
+      ``vlim`` at ± the maximum absolute value of the data and then chooses
+      contours within those bounds.
+
+    * When ``vlim`` but not a list of ``contours`` is passed, MNE chooses
+      contours to be within the ``vlim``.
+
+    * When a list of ``contours`` but not ``vlim`` is passed, MNE chooses
+      ``vlim`` to encompass the ``contours`` and the maximum absolute value of the
+      data.
+
+    * When both a list of ``contours`` and ``vlim`` are passed, MNE uses them
+      as-is.
+
     When ``time=="interactive"``, the figure will publish and subscribe to the
     following UI events:
 
@@ -2299,11 +2315,17 @@ def plot_evoked_topomap(
     _vlim = [
         _setup_vmin_vmax(data[:, i], *vlim, norm=merge_channels) for i in range(n_times)
     ]
-    _vlim = (np.min(_vlim), np.max(_vlim))
+    _vlim = [np.min(_vlim), np.max(_vlim)]
     cmap = _setup_cmap(cmap, n_axes=n_times, norm=_vlim[0] >= 0)
     # set up contours
     if not isinstance(contours, list | np.ndarray):
         _, contours = _set_contour_locator(*_vlim, contours)
+    else:
+        if vlim[0] is None and np.any(contours < _vlim[0]):
+            _vlim[0] = contours[0]
+        if vlim[1] is None and np.any(contours > _vlim[1]):
+            _vlim[1] = contours[-1]
+
     # prepare for main loop over times
     kwargs = dict(
         sensors=sensors,
@@ -3352,6 +3374,7 @@ def _set_contour_locator(vmin, vmax, contours):
         # correct number of bins is equal to contours + 1.
         locator = ticker.MaxNLocator(nbins=contours + 1)
         contours = locator.tick_values(vmin, vmax)
+        contours = contours[1:-1]
     return locator, contours