diff --git a/conftest.py b/conftest.py
index ddce5e0d593..da1dbc274fb 100644
--- a/conftest.py
+++ b/conftest.py
@@ -24,7 +24,7 @@ def pytest_runtest_setup(item):
 
 
 @pytest.fixture(autouse=True)
-def add_standard_imports(doctest_namespace):
+def add_standard_imports(doctest_namespace, tmpdir):
     import numpy as np
     import pandas as pd
 
@@ -33,3 +33,9 @@ def add_standard_imports(doctest_namespace):
     doctest_namespace["np"] = np
     doctest_namespace["pd"] = pd
     doctest_namespace["xr"] = xr
+
+    # always seed numpy.random to make the examples deterministic
+    np.random.seed(0)
+
+    # always switch to the temporary directory, so files get written there
+    tmpdir.chdir()
diff --git a/xarray/backends/api.py b/xarray/backends/api.py
index 74864478842..9f45474e7e7 100644
--- a/xarray/backends/api.py
+++ b/xarray/backends/api.py
@@ -1189,6 +1189,17 @@ def save_mfdataset(
 
     Save a dataset into one netCDF per year of data:
 
+    >>> ds = xr.Dataset(
+    ...     {"a": ("time", np.linspace(0, 1, 48))},
+    ...     coords={"time": pd.date_range("2010-01-01", freq="M", periods=48)},
+    ... )
+    >>> ds
+    <xarray.Dataset>
+    Dimensions:  (time: 48)
+    Coordinates:
+      * time     (time) datetime64[ns] 2010-01-31 2010-02-28 ... 2013-12-31
+    Data variables:
+        a        (time) float64 0.0 0.02128 0.04255 0.06383 ... 0.9574 0.9787 1.0
     >>> years, datasets = zip(*ds.groupby("time.year"))
     >>> paths = ["%s.nc" % y for y in years]
     >>> xr.save_mfdataset(datasets, paths)
diff --git a/xarray/coding/cftime_offsets.py b/xarray/coding/cftime_offsets.py
index 5ca4f5f6df3..99c7d041671 100644
--- a/xarray/coding/cftime_offsets.py
+++ b/xarray/coding/cftime_offsets.py
@@ -939,7 +939,7 @@ def cftime_range(
     >>> xr.cftime_range(start="2000", periods=6, freq="2MS", calendar="noleap")
     CFTimeIndex([2000-01-01 00:00:00, 2000-03-01 00:00:00, 2000-05-01 00:00:00,
                  2000-07-01 00:00:00, 2000-09-01 00:00:00, 2000-11-01 00:00:00],
-                dtype='object')
+                dtype='object', length=6, calendar='noleap')
 
     As in the standard pandas function, three of the ``start``, ``end``,
     ``periods``, or ``freq`` arguments must be specified at a given time, with
diff --git a/xarray/coding/cftimeindex.py b/xarray/coding/cftimeindex.py
index 85c6ee0809c..fe3f760f4a9 100644
--- a/xarray/coding/cftimeindex.py
+++ b/xarray/coding/cftimeindex.py
@@ -523,9 +523,11 @@ def shift(self, n, freq):
         --------
         >>> index = xr.cftime_range("2000", periods=1, freq="M")
         >>> index
-        CFTimeIndex([2000-01-31 00:00:00], dtype='object')
+        CFTimeIndex([2000-01-31 00:00:00],
+                    dtype='object', length=1, calendar='gregorian')
         >>> index.shift(1, "M")
-        CFTimeIndex([2000-02-29 00:00:00], dtype='object')
+        CFTimeIndex([2000-02-29 00:00:00],
+                    dtype='object', length=1, calendar='gregorian')
         """
         from .cftime_offsets import to_offset
 
@@ -611,7 +613,8 @@ def to_datetimeindex(self, unsafe=False):
         >>> import xarray as xr
         >>> times = xr.cftime_range("2000", periods=2, calendar="gregorian")
         >>> times
-        CFTimeIndex([2000-01-01 00:00:00, 2000-01-02 00:00:00], dtype='object')
+        CFTimeIndex([2000-01-01 00:00:00, 2000-01-02 00:00:00],
+                    dtype='object', length=2, calendar='gregorian')
         >>> times.to_datetimeindex()
         DatetimeIndex(['2000-01-01', '2000-01-02'], dtype='datetime64[ns]', freq=None)
         """
diff --git a/xarray/coding/strings.py b/xarray/coding/strings.py
index dfe0175947c..e16e983fd8a 100644
--- a/xarray/coding/strings.py
+++ b/xarray/coding/strings.py
@@ -199,9 +199,9 @@ class StackedBytesArray(indexing.ExplicitlyIndexedNDArrayMixin):
     """Wrapper around array-like objects to create a new indexable object where
     values, when accessed, are automatically stacked along the last dimension.
 
-    >>> StackedBytesArray(np.array(["a", "b", "c"]))[:]
-    array('abc',
-          dtype='|S3')
+    >>> indexer = indexing.BasicIndexer((slice(None),))
+    >>> StackedBytesArray(np.array(["a", "b", "c"], dtype="S1"))[indexer]
+    array(b'abc', dtype='|S3')
     """
 
     def __init__(self, array):
diff --git a/xarray/conventions.py b/xarray/conventions.py
index da69ce52527..bb0b92c77a1 100644
--- a/xarray/conventions.py
+++ b/xarray/conventions.py
@@ -24,10 +24,11 @@ class NativeEndiannessArray(indexing.ExplicitlyIndexedNDArrayMixin):
     >>> x.dtype
     dtype('>i2')
 
-    >>> NativeEndianArray(x).dtype
+    >>> NativeEndiannessArray(x).dtype
     dtype('int16')
 
-    >>> NativeEndianArray(x)[:].dtype
+    >>> indexer = indexing.BasicIndexer((slice(None),))
+    >>> NativeEndiannessArray(x)[indexer].dtype
     dtype('int16')
     """
 
@@ -53,12 +54,13 @@ class BoolTypeArray(indexing.ExplicitlyIndexedNDArrayMixin):
     >>> x = np.array([1, 0, 1, 1, 0], dtype="i1")
 
     >>> x.dtype
-    dtype('>i2')
+    dtype('int8')
 
     >>> BoolTypeArray(x).dtype
     dtype('bool')
 
-    >>> BoolTypeArray(x)[:].dtype
+    >>> indexer = indexing.BasicIndexer((slice(None),))
+    >>> BoolTypeArray(x)[indexer].dtype
     dtype('bool')
     """
 
diff --git a/xarray/core/accessor_dt.py b/xarray/core/accessor_dt.py
index a4ec7a2c30e..96ef0f3b5de 100644
--- a/xarray/core/accessor_dt.py
+++ b/xarray/core/accessor_dt.py
@@ -249,30 +249,30 @@ class DatetimeAccessor(Properties):
     >>> ts
     <xarray.DataArray (time: 10)>
     array(['2000-01-01T00:00:00.000000000', '2000-01-02T00:00:00.000000000',
-        '2000-01-03T00:00:00.000000000', '2000-01-04T00:00:00.000000000',
-        '2000-01-05T00:00:00.000000000', '2000-01-06T00:00:00.000000000',
-        '2000-01-07T00:00:00.000000000', '2000-01-08T00:00:00.000000000',
-        '2000-01-09T00:00:00.000000000', '2000-01-10T00:00:00.000000000'],
-        dtype='datetime64[ns]')
+           '2000-01-03T00:00:00.000000000', '2000-01-04T00:00:00.000000000',
+           '2000-01-05T00:00:00.000000000', '2000-01-06T00:00:00.000000000',
+           '2000-01-07T00:00:00.000000000', '2000-01-08T00:00:00.000000000',
+           '2000-01-09T00:00:00.000000000', '2000-01-10T00:00:00.000000000'],
+          dtype='datetime64[ns]')
     Coordinates:
-    * time     (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2000-01-10
-    >>> ts.dt
-    <xarray.core.accessor_dt.DatetimeAccessor object at 0x118b54d68>
+      * time     (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2000-01-10
+    >>> ts.dt  # doctest: +ELLIPSIS
+    <xarray.core.accessor_dt.DatetimeAccessor object at 0x...>
     >>> ts.dt.dayofyear
     <xarray.DataArray 'dayofyear' (time: 10)>
     array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
     Coordinates:
-    * time     (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2000-01-10
+      * time     (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2000-01-10
     >>> ts.dt.quarter
     <xarray.DataArray 'quarter' (time: 10)>
     array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
     Coordinates:
-    * time     (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2000-01-10
+      * time     (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2000-01-10
 
     """
 
     def strftime(self, date_format):
-        '''
+        """
         Return an array of formatted strings specified by date_format, which
         supports the same string format as the python standard library. Details
         of the string format can be found in `python string format doc
@@ -290,13 +290,12 @@ def strftime(self, date_format):
 
         Examples
         --------
+        >>> import datetime
         >>> rng = xr.Dataset({"time": datetime.datetime(2000, 1, 1)})
         >>> rng["time"].dt.strftime("%B %d, %Y, %r")
         <xarray.DataArray 'strftime' ()>
         array('January 01, 2000, 12:00:00 AM', dtype=object)
         """
-
-        '''
         obj_type = type(self._obj)
 
         result = _strftime(self._obj.data, date_format)
@@ -398,32 +397,32 @@ class TimedeltaAccessor(Properties):
     >>> ts
     <xarray.DataArray (time: 20)>
     array([ 86400000000000, 108000000000000, 129600000000000, 151200000000000,
-        172800000000000, 194400000000000, 216000000000000, 237600000000000,
-        259200000000000, 280800000000000, 302400000000000, 324000000000000,
-        345600000000000, 367200000000000, 388800000000000, 410400000000000,
-        432000000000000, 453600000000000, 475200000000000, 496800000000000],
-        dtype='timedelta64[ns]')
+           172800000000000, 194400000000000, 216000000000000, 237600000000000,
+           259200000000000, 280800000000000, 302400000000000, 324000000000000,
+           345600000000000, 367200000000000, 388800000000000, 410400000000000,
+           432000000000000, 453600000000000, 475200000000000, 496800000000000],
+          dtype='timedelta64[ns]')
     Coordinates:
-    * time     (time) timedelta64[ns] 1 days 00:00:00 ... 5 days 18:00:00
-    >>> ts.dt
-    <xarray.core.accessor_dt.TimedeltaAccessor object at 0x109a27d68>
+      * time     (time) timedelta64[ns] 1 days 00:00:00 ... 5 days 18:00:00
+    >>> ts.dt  # doctest: +ELLIPSIS
+    <xarray.core.accessor_dt.TimedeltaAccessor object at 0x...>
     >>> ts.dt.days
     <xarray.DataArray 'days' (time: 20)>
     array([1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5])
     Coordinates:
-    * time     (time) timedelta64[ns] 1 days 00:00:00 ... 5 days 18:00:00
+      * time     (time) timedelta64[ns] 1 days 00:00:00 ... 5 days 18:00:00
     >>> ts.dt.microseconds
     <xarray.DataArray 'microseconds' (time: 20)>
     array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
     Coordinates:
-    * time     (time) timedelta64[ns] 1 days 00:00:00 ... 5 days 18:00:00
+      * time     (time) timedelta64[ns] 1 days 00:00:00 ... 5 days 18:00:00
     >>> ts.dt.seconds
     <xarray.DataArray 'seconds' (time: 20)>
     array([    0, 21600, 43200, 64800,     0, 21600, 43200, 64800,     0,
-        21600, 43200, 64800,     0, 21600, 43200, 64800,     0, 21600,
-        43200, 64800])
+           21600, 43200, 64800,     0, 21600, 43200, 64800,     0, 21600,
+           43200, 64800])
     Coordinates:
-    * time     (time) timedelta64[ns] 1 days 00:00:00 ... 5 days 18:00:00
+      * time     (time) timedelta64[ns] 1 days 00:00:00 ... 5 days 18:00:00
     """
 
     days = Properties._tslib_field_accessor(
diff --git a/xarray/core/accessor_str.py b/xarray/core/accessor_str.py
index a845ce1f642..faf212e86ee 100644
--- a/xarray/core/accessor_str.py
+++ b/xarray/core/accessor_str.py
@@ -68,7 +68,7 @@ class StringAccessor:
     for applicable DataArrays.
 
         >>> da = xr.DataArray(["some", "text", "in", "an", "array"])
-        >>> ds.str.len()
+        >>> da.str.len()
         <xarray.DataArray (dim_0: 5)>
         array([4, 4, 2, 2, 5])
         Dimensions without coordinates: dim_0
diff --git a/xarray/core/alignment.py b/xarray/core/alignment.py
index 23a3cc719a8..319c8a9a7a0 100644
--- a/xarray/core/alignment.py
+++ b/xarray/core/alignment.py
@@ -140,30 +140,30 @@ def align(
     array([[25, 35],
            [10, 24]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0
+      * lon      (lon) float64 100.0 120.0
 
     >>> y
     <xarray.DataArray (lat: 2, lon: 2)>
     array([[20,  5],
            [ 7, 13]])
     Coordinates:
-    * lat      (lat) float64 35.0 42.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 42.0
+      * lon      (lon) float64 100.0 120.0
 
     >>> a, b = xr.align(x, y)
     >>> a
     <xarray.DataArray (lat: 1, lon: 2)>
     array([[25, 35]])
     Coordinates:
-    * lat      (lat) float64 35.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0
+      * lon      (lon) float64 100.0 120.0
     >>> b
     <xarray.DataArray (lat: 1, lon: 2)>
     array([[20,  5]])
     Coordinates:
-    * lat      (lat) float64 35.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0
+      * lon      (lon) float64 100.0 120.0
 
     >>> a, b = xr.align(x, y, join="outer")
     >>> a
@@ -172,16 +172,16 @@ def align(
            [10., 24.],
            [nan, nan]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0
     >>> b
     <xarray.DataArray (lat: 3, lon: 2)>
     array([[20.,  5.],
            [nan, nan],
            [ 7., 13.]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0
 
     >>> a, b = xr.align(x, y, join="outer", fill_value=-999)
     >>> a
@@ -190,16 +190,16 @@ def align(
            [  10,   24],
            [-999, -999]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0
     >>> b
     <xarray.DataArray (lat: 3, lon: 2)>
     array([[  20,    5],
            [-999, -999],
            [   7,   13]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0
 
     >>> a, b = xr.align(x, y, join="left")
     >>> a
@@ -207,15 +207,15 @@ def align(
     array([[25, 35],
            [10, 24]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0
+      * lon      (lon) float64 100.0 120.0
     >>> b
     <xarray.DataArray (lat: 2, lon: 2)>
     array([[20.,  5.],
            [nan, nan]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0
+      * lon      (lon) float64 100.0 120.0
 
     >>> a, b = xr.align(x, y, join="right")
     >>> a
@@ -223,15 +223,15 @@ def align(
     array([[25., 35.],
            [nan, nan]])
     Coordinates:
-    * lat      (lat) float64 35.0 42.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 42.0
+      * lon      (lon) float64 100.0 120.0
     >>> b
     <xarray.DataArray (lat: 2, lon: 2)>
     array([[20,  5],
            [ 7, 13]])
     Coordinates:
-    * lat      (lat) float64 35.0 42.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 42.0
+      * lon      (lon) float64 100.0 120.0
 
     >>> a, b = xr.align(x, y, join="exact")
     Traceback (most recent call last):
@@ -245,15 +245,15 @@ def align(
     array([[25, 35],
            [10, 24]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0
+      * lon      (lon) float64 100.0 120.0
     >>> b
     <xarray.DataArray (lat: 2, lon: 2)>
     array([[20,  5],
            [ 7, 13]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0
+      * lon      (lon) float64 100.0 120.0
 
     """
     if indexes is None:
@@ -693,30 +693,24 @@ def broadcast(*args, exclude=None):
     >>> a
     <xarray.DataArray (x: 3)>
     array([1, 2, 3])
-    Coordinates:
-      * x        (x) int64 0 1 2
+    Dimensions without coordinates: x
     >>> b
     <xarray.DataArray (y: 2)>
     array([5, 6])
-    Coordinates:
-      * y        (y) int64 0 1
+    Dimensions without coordinates: y
     >>> a2, b2 = xr.broadcast(a, b)
     >>> a2
     <xarray.DataArray (x: 3, y: 2)>
     array([[1, 1],
            [2, 2],
            [3, 3]])
-    Coordinates:
-      * x        (x) int64 0 1 2
-      * y        (y) int64 0 1
+    Dimensions without coordinates: x, y
     >>> b2
     <xarray.DataArray (x: 3, y: 2)>
     array([[5, 6],
            [5, 6],
            [5, 6]])
-    Coordinates:
-      * y        (y) int64 0 1
-      * x        (x) int64 0 1 2
+    Dimensions without coordinates: x, y
 
     Fill out the dimensions of all data variables in a dataset:
 
@@ -725,9 +719,7 @@ def broadcast(*args, exclude=None):
     >>> ds2
     <xarray.Dataset>
     Dimensions:  (x: 3, y: 2)
-    Coordinates:
-      * x        (x) int64 0 1 2
-      * y        (y) int64 0 1
+    Dimensions without coordinates: x, y
     Data variables:
         a        (x, y) int64 1 1 2 2 3 3
         b        (x, y) int64 5 6 5 6 5 6
diff --git a/xarray/core/combine.py b/xarray/core/combine.py
index 7fc9e101cd2..d9ce3def673 100644
--- a/xarray/core/combine.py
+++ b/xarray/core/combine.py
@@ -433,22 +433,48 @@ def combine_nested(
     into 4 parts, 2 each along both the x and y axes, requires organising the
     datasets into a doubly-nested list, e.g:
 
+    >>> x1y1 = xr.Dataset(
+    ...     {
+    ...         "temperature": (("x", "y"), np.random.randn(2, 2)),
+    ...         "precipitation": (("x", "y"), np.random.randn(2, 2)),
+    ...     }
+    ... )
     >>> x1y1
     <xarray.Dataset>
-    Dimensions:         (x: 2, y: 2)
+    Dimensions:        (x: 2, y: 2)
     Dimensions without coordinates: x, y
     Data variables:
-      temperature       (x, y) float64 11.04 23.57 20.77 ...
-      precipitation     (x, y) float64 5.904 2.453 3.404 ...
+        temperature    (x, y) float64 1.764 0.4002 0.9787 2.241
+        precipitation  (x, y) float64 1.868 -0.9773 0.9501 -0.1514
+    >>> x1y2 = xr.Dataset(
+    ...     {
+    ...         "temperature": (("x", "y"), np.random.randn(2, 2)),
+    ...         "precipitation": (("x", "y"), np.random.randn(2, 2)),
+    ...     }
+    ... )
+    >>> x2y1 = xr.Dataset(
+    ...     {
+    ...         "temperature": (("x", "y"), np.random.randn(2, 2)),
+    ...         "precipitation": (("x", "y"), np.random.randn(2, 2)),
+    ...     }
+    ... )
+    >>> x2y2 = xr.Dataset(
+    ...     {
+    ...         "temperature": (("x", "y"), np.random.randn(2, 2)),
+    ...         "precipitation": (("x", "y"), np.random.randn(2, 2)),
+    ...     }
+    ... )
+
 
     >>> ds_grid = [[x1y1, x1y2], [x2y1, x2y2]]
     >>> combined = xr.combine_nested(ds_grid, concat_dim=["x", "y"])
+    >>> combined
     <xarray.Dataset>
-    Dimensions:         (x: 4, y: 4)
+    Dimensions:        (x: 4, y: 4)
     Dimensions without coordinates: x, y
     Data variables:
-      temperature       (x, y) float64 11.04 23.57 20.77 ...
-      precipitation     (x, y) float64 5.904 2.453 3.404 ...
+        temperature    (x, y) float64 1.764 0.4002 -0.1032 ... 0.04576 -0.1872
+        precipitation  (x, y) float64 1.868 -0.9773 0.761 ... -0.7422 0.1549 0.3782
 
     ``manual_combine`` can also be used to explicitly merge datasets with
     different variables. For example if we have 4 datasets, which are divided
@@ -456,28 +482,35 @@ def combine_nested(
     to ``concat_dim`` to specify the dimension of the nested list over which
     we wish to use ``merge`` instead of ``concat``:
 
+    >>> t1temp = xr.Dataset({"temperature": ("t", np.random.randn(5))})
     >>> t1temp
     <xarray.Dataset>
-    Dimensions:         (t: 5)
+    Dimensions:      (t: 5)
     Dimensions without coordinates: t
     Data variables:
-      temperature       (t) float64 11.04 23.57 20.77 ...
+        temperature  (t) float64 -0.8878 -1.981 -0.3479 0.1563 1.23
 
+    >>> t1precip = xr.Dataset({"precipitation": ("t", np.random.randn(5))})
     >>> t1precip
     <xarray.Dataset>
-    Dimensions:         (t: 5)
+    Dimensions:        (t: 5)
     Dimensions without coordinates: t
     Data variables:
-      precipitation     (t) float64 5.904 2.453 3.404 ...
+        precipitation  (t) float64 1.202 -0.3873 -0.3023 -1.049 -1.42
+
+    >>> t2temp = xr.Dataset({"temperature": ("t", np.random.randn(5))})
+    >>> t2precip = xr.Dataset({"precipitation": ("t", np.random.randn(5))})
+
 
     >>> ds_grid = [[t1temp, t1precip], [t2temp, t2precip]]
     >>> combined = xr.combine_nested(ds_grid, concat_dim=["t", None])
+    >>> combined
     <xarray.Dataset>
-    Dimensions:         (t: 10)
+    Dimensions:        (t: 10)
     Dimensions without coordinates: t
     Data variables:
-      temperature       (t) float64 11.04 23.57 20.77 ...
-      precipitation     (t) float64 5.904 2.453 3.404 ...
+        temperature    (t) float64 -0.8878 -1.981 -0.3479 ... -0.5097 -0.4381 -1.253
+        precipitation  (t) float64 1.202 -0.3873 -0.3023 ... -0.2127 -0.8955 0.3869
 
     See also
     --------
@@ -648,71 +681,71 @@ def combine_by_coords(
     <xarray.Dataset>
     Dimensions:        (x: 3, y: 2)
     Coordinates:
-    * y              (y) int64 0 1
-    * x              (x) int64 10 20 30
+      * y              (y) int64 0 1
+      * x              (x) int64 10 20 30
     Data variables:
-        temperature    (y, x) float64 1.654 10.63 7.015 2.543 13.93 9.436
-        precipitation  (y, x) float64 0.2136 0.9974 0.7603 0.4679 0.3115 0.945
+        temperature    (y, x) float64 10.98 14.3 12.06 10.9 8.473 12.92
+        precipitation  (y, x) float64 0.4376 0.8918 0.9637 0.3834 0.7917 0.5289
 
     >>> x2
     <xarray.Dataset>
     Dimensions:        (x: 3, y: 2)
     Coordinates:
-    * y              (y) int64 2 3
-    * x              (x) int64 10 20 30
+      * y              (y) int64 2 3
+      * x              (x) int64 10 20 30
     Data variables:
-        temperature    (y, x) float64 9.341 0.1251 6.269 7.709 8.82 2.316
-        precipitation  (y, x) float64 0.1728 0.1178 0.03018 0.6509 0.06938 0.3792
+        temperature    (y, x) float64 11.36 18.51 1.421 1.743 0.4044 16.65
+        precipitation  (y, x) float64 0.7782 0.87 0.9786 0.7992 0.4615 0.7805
 
     >>> x3
     <xarray.Dataset>
     Dimensions:        (x: 3, y: 2)
     Coordinates:
-    * y              (y) int64 2 3
-    * x              (x) int64 40 50 60
+      * y              (y) int64 2 3
+      * x              (x) int64 40 50 60
     Data variables:
-        temperature    (y, x) float64 2.789 2.446 6.551 12.46 2.22 15.96
-        precipitation  (y, x) float64 0.4804 0.1902 0.2457 0.6125 0.4654 0.5953
+        temperature    (y, x) float64 2.365 12.8 2.867 18.89 10.44 8.293
+        precipitation  (y, x) float64 0.2646 0.7742 0.4562 0.5684 0.01879 0.6176
 
     >>> xr.combine_by_coords([x2, x1])
     <xarray.Dataset>
     Dimensions:        (x: 3, y: 4)
     Coordinates:
-    * x              (x) int64 10 20 30
-    * y              (y) int64 0 1 2 3
+      * x              (x) int64 10 20 30
+      * y              (y) int64 0 1 2 3
     Data variables:
-        temperature    (y, x) float64 1.654 10.63 7.015 2.543 ... 7.709 8.82 2.316
-        precipitation  (y, x) float64 0.2136 0.9974 0.7603 ... 0.6509 0.06938 0.3792
+        temperature    (y, x) float64 10.98 14.3 12.06 10.9 ... 1.743 0.4044 16.65
+        precipitation  (y, x) float64 0.4376 0.8918 0.9637 ... 0.7992 0.4615 0.7805
 
     >>> xr.combine_by_coords([x3, x1])
     <xarray.Dataset>
     Dimensions:        (x: 6, y: 4)
     Coordinates:
-    * x              (x) int64 10 20 30 40 50 60
-    * y              (y) int64 0 1 2 3
+      * x              (x) int64 10 20 30 40 50 60
+      * y              (y) int64 0 1 2 3
     Data variables:
-        temperature    (y, x) float64 1.654 10.63 7.015 nan ... nan 12.46 2.22 15.96
-        precipitation  (y, x) float64 0.2136 0.9974 0.7603 ... 0.6125 0.4654 0.5953
+        temperature    (y, x) float64 10.98 14.3 12.06 nan ... nan 18.89 10.44 8.293
+        precipitation  (y, x) float64 0.4376 0.8918 0.9637 ... 0.5684 0.01879 0.6176
 
     >>> xr.combine_by_coords([x3, x1], join="override")
     <xarray.Dataset>
     Dimensions:        (x: 3, y: 4)
     Coordinates:
-    * x              (x) int64 10 20 30
-    * y              (y) int64 0 1 2 3
+      * x              (x) int64 10 20 30
+      * y              (y) int64 0 1 2 3
     Data variables:
-    temperature    (y, x) float64 1.654 10.63 7.015 2.543 ... 12.46 2.22 15.96
-    precipitation  (y, x) float64 0.2136 0.9974 0.7603 ... 0.6125 0.4654 0.5953
+        temperature    (y, x) float64 10.98 14.3 12.06 10.9 ... 18.89 10.44 8.293
+        precipitation  (y, x) float64 0.4376 0.8918 0.9637 ... 0.5684 0.01879 0.6176
 
     >>> xr.combine_by_coords([x1, x2, x3])
     <xarray.Dataset>
     Dimensions:        (x: 6, y: 4)
     Coordinates:
-    * x              (x) int64 10 20 30 40 50 60
-    * y              (y) int64 0 1 2 3
+      * x              (x) int64 10 20 30 40 50 60
+      * y              (y) int64 0 1 2 3
     Data variables:
-    temperature    (y, x) float64 1.654 10.63 7.015 nan ... 12.46 2.22 15.96
-    precipitation  (y, x) float64 0.2136 0.9974 0.7603 ... 0.6125 0.4654 0.5953
+        temperature    (y, x) float64 10.98 14.3 12.06 nan ... 18.89 10.44 8.293
+        precipitation  (y, x) float64 0.4376 0.8918 0.9637 ... 0.5684 0.01879 0.6176
     """
 
     # Group by data vars
diff --git a/xarray/core/common.py b/xarray/core/common.py
index 38803f821d4..b48a2f56a0d 100644
--- a/xarray/core/common.py
+++ b/xarray/core/common.py
@@ -422,12 +422,12 @@ def assign_coords(self, coords=None, **coords_kwargs):
         ... )
         >>> da
         <xarray.DataArray (lon: 4)>
-        array([0.28298 , 0.667347, 0.657938, 0.177683])
+        array([0.5488135 , 0.71518937, 0.60276338, 0.54488318])
         Coordinates:
           * lon      (lon) int64 358 359 0 1
         >>> da.assign_coords(lon=(((da.lon + 180) % 360) - 180))
         <xarray.DataArray (lon: 4)>
-        array([0.28298 , 0.667347, 0.657938, 0.177683])
+        array([0.5488135 , 0.71518937, 0.60276338, 0.54488318])
         Coordinates:
           * lon      (lon) int64 -2 -1 0 1
 
@@ -435,7 +435,7 @@ def assign_coords(self, coords=None, **coords_kwargs):
 
         >>> da.assign_coords({"lon": (((da.lon + 180) % 360) - 180)})
         <xarray.DataArray (lon: 4)>
-        array([0.28298 , 0.667347, 0.657938, 0.177683])
+        array([0.5488135 , 0.71518937, 0.60276338, 0.54488318])
         Coordinates:
           * lon      (lon) int64 -2 -1 0 1
 
@@ -444,7 +444,7 @@ def assign_coords(self, coords=None, **coords_kwargs):
         >>> lon_2 = np.array([300, 289, 0, 1])
         >>> da.assign_coords(lon_2=("lon", lon_2))
         <xarray.DataArray (lon: 4)>
-        array([0.28298 , 0.667347, 0.657938, 0.177683])
+        array([0.5488135 , 0.71518937, 0.60276338, 0.54488318])
         Coordinates:
           * lon      (lon) int64 358 359 0 1
             lon_2    (lon) int64 300 289 0 1
@@ -532,17 +532,23 @@ def pipe(
         Use ``.pipe`` when chaining together functions that expect
         xarray or pandas objects, e.g., instead of writing
 
-        >>> f(g(h(ds), arg1=a), arg2=b, arg3=c)
+        .. code:: python
+
+            f(g(h(ds), arg1=a), arg2=b, arg3=c)
 
         You can write
 
-        >>> (ds.pipe(h).pipe(g, arg1=a).pipe(f, arg2=b, arg3=c))
+        .. code:: python
+
+            (ds.pipe(h).pipe(g, arg1=a).pipe(f, arg2=b, arg3=c))
 
         If you have a function that takes the data as (say) the second
         argument, pass a tuple indicating which keyword expects the
         data. For example, suppose ``f`` takes its data as ``arg2``:
 
-        >>> (ds.pipe(h).pipe(g, arg1=a).pipe((f, "arg2"), arg1=a, arg3=c))
+        .. code:: python
+
+            (ds.pipe(h).pipe(g, arg1=a).pipe((f, "arg2"), arg1=a, arg3=c))
 
         Examples
         --------
@@ -563,11 +569,11 @@ def pipe(
         <xarray.Dataset>
         Dimensions:        (lat: 2, lon: 2)
         Coordinates:
-        * lat            (lat) int64 10 20
-        * lon            (lon) int64 150 160
+          * lat            (lat) int64 10 20
+          * lon            (lon) int64 150 160
         Data variables:
-            temperature_c  (lat, lon) float64 14.53 11.85 19.27 16.37
-            precipitation  (lat, lon) float64 0.7315 0.7189 0.8481 0.4671
+            temperature_c  (lat, lon) float64 10.98 14.3 12.06 10.9
+            precipitation  (lat, lon) float64 0.4237 0.6459 0.4376 0.8918
 
         >>> def adder(data, arg):
         ...     return data + arg
@@ -582,21 +588,21 @@ def pipe(
         <xarray.Dataset>
         Dimensions:        (lat: 2, lon: 2)
         Coordinates:
-        * lon            (lon) int64 150 160
-        * lat            (lat) int64 10 20
+          * lat            (lat) int64 10 20
+          * lon            (lon) int64 150 160
         Data variables:
-            temperature_c  (lat, lon) float64 16.53 13.85 21.27 18.37
-            precipitation  (lat, lon) float64 2.731 2.719 2.848 2.467
+            temperature_c  (lat, lon) float64 12.98 16.3 14.06 12.9
+            precipitation  (lat, lon) float64 2.424 2.646 2.438 2.892
 
         >>> x.pipe(adder, arg=2)
         <xarray.Dataset>
         Dimensions:        (lat: 2, lon: 2)
         Coordinates:
-        * lon            (lon) int64 150 160
-        * lat            (lat) int64 10 20
+          * lat            (lat) int64 10 20
+          * lon            (lon) int64 150 160
         Data variables:
-            temperature_c  (lat, lon) float64 16.53 13.85 21.27 18.37
-            precipitation  (lat, lon) float64 2.731 2.719 2.848 2.467
+            temperature_c  (lat, lon) float64 12.98 16.3 14.06 12.9
+            precipitation  (lat, lon) float64 2.424 2.646 2.438 2.892
 
         >>> (
         ...     x.pipe(adder, arg=2)
@@ -606,11 +612,11 @@ def pipe(
         <xarray.Dataset>
         Dimensions:        (lat: 2, lon: 2)
         Coordinates:
-        * lon            (lon) int64 150 160
-        * lat            (lat) int64 10 20
+          * lat            (lat) int64 10 20
+          * lon            (lon) int64 150 160
         Data variables:
-            temperature_c  (lat, lon) float64 14.53 11.85 19.27 16.37
-            precipitation  (lat, lon) float64 0.7315 0.7189 0.8481 0.4671
+            temperature_c  (lat, lon) float64 10.98 14.3 12.06 10.9
+            precipitation  (lat, lon) float64 0.4237 0.6459 0.4376 0.8918
 
         See Also
         --------
@@ -660,15 +666,16 @@ def groupby(self, group, squeeze: bool = True, restore_coord_dims: bool = None):
         ... )
         >>> da
         <xarray.DataArray (time: 1827)>
-        array([0.000e+00, 1.000e+00, 2.000e+00, ..., 1.824e+03, 1.825e+03, 1.826e+03])
+        array([0.000e+00, 1.000e+00, 2.000e+00, ..., 1.824e+03, 1.825e+03,
+               1.826e+03])
         Coordinates:
-          * time     (time) datetime64[ns] 2000-01-01 2000-01-02 2000-01-03 ...
+          * time     (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2004-12-31
         >>> da.groupby("time.dayofyear") - da.groupby("time.dayofyear").mean("time")
         <xarray.DataArray (time: 1827)>
         array([-730.8, -730.8, -730.8, ...,  730.2,  730.2,  730.5])
         Coordinates:
-          * time       (time) datetime64[ns] 2000-01-01 2000-01-02 2000-01-03 ...
-            dayofyear  (time) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ...
+          * time       (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2004-12-31
+            dayofyear  (time) int64 1 2 3 4 5 6 7 8 ... 359 360 361 362 363 364 365 366
 
         See Also
         --------
@@ -834,14 +841,14 @@ def rolling(
         ... )
         >>> da
         <xarray.DataArray (time: 12)>
-        array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7., 8.,   9.,  10.,  11.])
+        array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11.])
         Coordinates:
-          * time     (time) datetime64[ns] 1999-12-15 2000-01-15 2000-02-15 ...
+          * time     (time) datetime64[ns] 1999-12-15 2000-01-15 ... 2000-11-15
         >>> da.rolling(time=3, center=True).mean()
         <xarray.DataArray (time: 12)>
         array([nan,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., nan])
         Coordinates:
-          * time     (time) datetime64[ns] 1999-12-15 2000-01-15 2000-02-15 ...
+          * time     (time) datetime64[ns] 1999-12-15 2000-01-15 ... 2000-11-15
 
         Remove the NaNs using ``dropna()``:
 
@@ -849,7 +856,7 @@ def rolling(
         <xarray.DataArray (time: 10)>
         array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.])
         Coordinates:
-          * time     (time) datetime64[ns] 2000-01-15 2000-02-15 2000-03-15 ...
+          * time     (time) datetime64[ns] 2000-01-15 2000-02-15 ... 2000-10-15
 
         See Also
         --------
@@ -940,17 +947,24 @@ def coarsen(
         ...     dims="time",
         ...     coords={"time": pd.date_range("15/12/1999", periods=364)},
         ... )
-        >>> da
+        >>> da  # +doctest: ELLIPSIS
         <xarray.DataArray (time: 364)>
-        array([  0.      ,   1.002755,   2.00551 , ..., 361.99449 , 362.997245,
-               364.      ])
+        array([  0.        ,   1.00275482,   2.00550964,   3.00826446,
+                 4.01101928,   5.0137741 ,   6.01652893,   7.01928375,
+                 8.02203857,   9.02479339,  10.02754821,  11.03030303,
+        ...
+               356.98071625, 357.98347107, 358.9862259 , 359.98898072,
+               360.99173554, 361.99449036, 362.99724518, 364.        ])
         Coordinates:
           * time     (time) datetime64[ns] 1999-12-15 1999-12-16 ... 2000-12-12
-        >>>
-        >>> da.coarsen(time=3, boundary="trim").mean()
+        >>> da.coarsen(time=3, boundary="trim").mean()  # +doctest: ELLIPSIS
         <xarray.DataArray (time: 121)>
-        array([  1.002755,   4.011019,   7.019284,  ...,  358.986226,
-               361.99449 ])
+        array([  1.00275482,   4.01101928,   7.01928375,  10.02754821,
+                13.03581267,  16.04407713,  19.0523416 ,  22.06060606,
+                25.06887052,  28.07713499,  31.08539945,  34.09366391,
+        ...
+               349.96143251, 352.96969697, 355.97796143, 358.9862259 ,
+               361.99449036])
         Coordinates:
           * time     (time) datetime64[ns] 1999-12-16 1999-12-19 ... 2000-12-10
         >>>
@@ -1043,9 +1057,9 @@ def resample(
         ... )
         >>> da
         <xarray.DataArray (time: 12)>
-        array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7., 8.,   9.,  10.,  11.])
+        array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11.])
         Coordinates:
-          * time     (time) datetime64[ns] 1999-12-15 2000-01-15 2000-02-15 ...
+          * time     (time) datetime64[ns] 1999-12-15 2000-01-15 ... 2000-11-15
         >>> da.resample(time="QS-DEC").mean()
         <xarray.DataArray (time: 4)>
         array([ 1.,  4.,  7., 10.])
@@ -1054,11 +1068,16 @@ def resample(
 
         Upsample monthly time-series data to daily data:
 
-        >>> da.resample(time="1D").interpolate("linear")
+        >>> da.resample(time="1D").interpolate("linear")  # +doctest: ELLIPSIS
         <xarray.DataArray (time: 337)>
-        array([ 0.      ,  0.032258,  0.064516, ..., 10.935484, 10.967742, 11.      ])
+        array([ 0.        ,  0.03225806,  0.06451613,  0.09677419,  0.12903226,
+                0.16129032,  0.19354839,  0.22580645,  0.25806452,  0.29032258,
+                0.32258065,  0.35483871,  0.38709677,  0.41935484,  0.4516129 ,
+        ...
+               10.80645161, 10.83870968, 10.87096774, 10.90322581, 10.93548387,
+               10.96774194, 11.        ])
         Coordinates:
-          * time     (time) datetime64[ns] 1999-12-15 1999-12-16 1999-12-17 ...
+          * time     (time) datetime64[ns] 1999-12-15 1999-12-16 ... 2000-11-15
 
         Limit scope of upsampling method
 
@@ -1160,19 +1179,19 @@ def where(self, cond, other=dtypes.NA, drop: bool = False):
         >>> a
         <xarray.DataArray (x: 5, y: 5)>
         array([[ 0,  1,  2,  3,  4],
-            [ 5,  6,  7,  8,  9],
-            [10, 11, 12, 13, 14],
-            [15, 16, 17, 18, 19],
-            [20, 21, 22, 23, 24]])
+               [ 5,  6,  7,  8,  9],
+               [10, 11, 12, 13, 14],
+               [15, 16, 17, 18, 19],
+               [20, 21, 22, 23, 24]])
         Dimensions without coordinates: x, y
 
         >>> a.where(a.x + a.y < 4)
         <xarray.DataArray (x: 5, y: 5)>
-        array([[  0.,   1.,   2.,   3.,  nan],
-               [  5.,   6.,   7.,  nan,  nan],
-               [ 10.,  11.,  nan,  nan,  nan],
-               [ 15.,  nan,  nan,  nan,  nan],
-               [ nan,  nan,  nan,  nan,  nan]])
+        array([[ 0.,  1.,  2.,  3., nan],
+               [ 5.,  6.,  7., nan, nan],
+               [10., 11., nan, nan, nan],
+               [15., nan, nan, nan, nan],
+               [nan, nan, nan, nan, nan]])
         Dimensions without coordinates: x, y
 
         >>> a.where(a.x + a.y < 5, -1)
@@ -1186,18 +1205,18 @@ def where(self, cond, other=dtypes.NA, drop: bool = False):
 
         >>> a.where(a.x + a.y < 4, drop=True)
         <xarray.DataArray (x: 4, y: 4)>
-        array([[  0.,   1.,   2.,   3.],
-               [  5.,   6.,   7.,  nan],
-               [ 10.,  11.,  nan,  nan],
-               [ 15.,  nan,  nan,  nan]])
+        array([[ 0.,  1.,  2.,  3.],
+               [ 5.,  6.,  7., nan],
+               [10., 11., nan, nan],
+               [15., nan, nan, nan]])
         Dimensions without coordinates: x, y
 
         >>> a.where(lambda x: x.x + x.y < 4, drop=True)
         <xarray.DataArray (x: 4, y: 4)>
-        array([[  0.,   1.,   2.,   3.],
-               [  5.,   6.,   7.,  nan],
-               [ 10.,  11.,  nan,  nan],
-               [ 15.,  nan,  nan,  nan]])
+        array([[ 0.,  1.,  2.,  3.],
+               [ 5.,  6.,  7., nan],
+               [10., 11., nan, nan],
+               [15., nan, nan, nan]])
         Dimensions without coordinates: x, y
 
         See also
@@ -1393,40 +1412,40 @@ def full_like(other, fill_value, dtype: DTypeLike = None):
     array([[0, 1, 2],
            [3, 4, 5]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     >>> xr.full_like(x, 1)
     <xarray.DataArray (lat: 2, lon: 3)>
     array([[1, 1, 1],
            [1, 1, 1]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     >>> xr.full_like(x, 0.5)
     <xarray.DataArray (lat: 2, lon: 3)>
     array([[0, 0, 0],
            [0, 0, 0]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     >>> xr.full_like(x, 0.5, dtype=np.double)
     <xarray.DataArray (lat: 2, lon: 3)>
     array([[0.5, 0.5, 0.5],
            [0.5, 0.5, 0.5]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     >>> xr.full_like(x, np.nan, dtype=np.double)
     <xarray.DataArray (lat: 2, lon: 3)>
     array([[nan, nan, nan],
            [nan, nan, nan]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     >>> ds = xr.Dataset(
     ...     {"a": ("x", [3, 5, 2]), "b": ("x", [9, 1, 0])}, coords={"x": [2, 4, 6]}
@@ -1552,24 +1571,24 @@ def zeros_like(other, dtype: DTypeLike = None):
     array([[0, 1, 2],
            [3, 4, 5]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     >>> xr.zeros_like(x)
     <xarray.DataArray (lat: 2, lon: 3)>
     array([[0, 0, 0],
            [0, 0, 0]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     >>> xr.zeros_like(x, dtype=float)
     <xarray.DataArray (lat: 2, lon: 3)>
     array([[0., 0., 0.],
            [0., 0., 0.]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     See also
     --------
@@ -1612,16 +1631,16 @@ def ones_like(other, dtype: DTypeLike = None):
     array([[0, 1, 2],
            [3, 4, 5]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     >>> xr.ones_like(x)
     <xarray.DataArray (lat: 2, lon: 3)>
     array([[1, 1, 1],
            [1, 1, 1]])
     Coordinates:
-    * lat      (lat) int64 1 2
-    * lon      (lon) int64 0 1 2
+      * lat      (lat) int64 1 2
+      * lon      (lon) int64 0 1 2
 
     See also
     --------
diff --git a/xarray/core/computation.py b/xarray/core/computation.py
index 507f12fe55e..85b88c4a976 100644
--- a/xarray/core/computation.py
+++ b/xarray/core/computation.py
@@ -894,7 +894,7 @@ def apply_ufunc(
     >>> array = xr.DataArray([1, 2, 3], coords=[("x", [0.1, 0.2, 0.3])])
     >>> magnitude(array, -array)
     <xarray.DataArray (x: 3)>
-    array([1.414214, 2.828427, 4.242641])
+    array([1.41421356, 2.82842712, 4.24264069])
     Coordinates:
       * x        (x) float64 0.1 0.2 0.3
 
@@ -1124,6 +1124,7 @@ def cov(da_a, da_b, dim=None, ddof=1):
 
     Examples
     --------
+    >>> from xarray import DataArray
     >>> da_a = DataArray(
     ...     np.array([[1, 2, 3], [0.1, 0.2, 0.3], [3.2, 0.6, 1.8]]),
     ...     dims=("space", "time"),
@@ -1161,7 +1162,7 @@ def cov(da_a, da_b, dim=None, ddof=1):
     array(-3.53055556)
     >>> xr.cov(da_a, da_b, dim="time")
     <xarray.DataArray (space: 3)>
-    array([ 0.2, -0.5,  1.69333333])
+    array([ 0.2       , -0.5       ,  1.69333333])
     Coordinates:
       * space    (space) <U2 'IA' 'IL' 'IN'
     """
@@ -1201,6 +1202,7 @@ def corr(da_a, da_b, dim=None):
 
     Examples
     --------
+    >>> from xarray import DataArray
     >>> da_a = DataArray(
     ...     np.array([[1, 2, 3], [0.1, 0.2, 0.3], [3.2, 0.6, 1.8]]),
     ...     dims=("space", "time"),
@@ -1332,8 +1334,10 @@ def dot(*arrays, dims=None, **kwargs):
     <xarray.DataArray (a: 3, b: 2, c: 2)>
     array([[[ 0,  1],
             [ 2,  3]],
+    <BLANKLINE>
            [[ 4,  5],
             [ 6,  7]],
+    <BLANKLINE>
            [[ 8,  9],
             [10, 11]]])
     Dimensions without coordinates: a, b, c
@@ -1477,13 +1481,13 @@ def where(cond, x, y):
     <xarray.DataArray 'sst' (lat: 10)>
     array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
     Coordinates:
-    * lat      (lat) int64 0 1 2 3 4 5 6 7 8 9
+      * lat      (lat) int64 0 1 2 3 4 5 6 7 8 9
 
     >>> xr.where(x < 0.5, x, x * 100)
     <xarray.DataArray 'sst' (lat: 10)>
     array([ 0. ,  0.1,  0.2,  0.3,  0.4, 50. , 60. , 70. , 80. , 90. ])
     Coordinates:
-    * lat      (lat) int64 0 1 2 3 4 5 6 7 8 9
+      * lat      (lat) int64 0 1 2 3 4 5 6 7 8 9
 
     >>> y = xr.DataArray(
     ...     0.1 * np.arange(9).reshape(3, 3),
@@ -1497,8 +1501,8 @@ def where(cond, x, y):
            [0.3, 0.4, 0.5],
            [0.6, 0.7, 0.8]])
     Coordinates:
-    * lat      (lat) int64 0 1 2
-    * lon      (lon) int64 10 11 12
+      * lat      (lat) int64 0 1 2
+      * lon      (lon) int64 10 11 12
 
     >>> xr.where(y.lat < 1, y, -1)
     <xarray.DataArray (lat: 3, lon: 3)>
@@ -1506,8 +1510,8 @@ def where(cond, x, y):
            [-1. , -1. , -1. ],
            [-1. , -1. , -1. ]])
     Coordinates:
-    * lat      (lat) int64 0 1 2
-    * lon      (lon) int64 10 11 12
+      * lat      (lat) int64 0 1 2
+      * lon      (lon) int64 10 11 12
 
     >>> cond = xr.DataArray([True, False], dims=["x"])
     >>> x = xr.DataArray([1, 2], dims=["y"])
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index 7b3af0f8391..86cb7ad988e 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -888,19 +888,19 @@ def copy(self, deep: bool = True, data: Any = None) -> "DataArray":
         <xarray.DataArray (x: 3)>
         array([1, 2, 3])
         Coordinates:
-        * x        (x) <U1 'a' 'b' 'c'
+          * x        (x) <U1 'a' 'b' 'c'
         >>> array_0 = array.copy(deep=False)
         >>> array_0[0] = 7
         >>> array_0
         <xarray.DataArray (x: 3)>
         array([7, 2, 3])
         Coordinates:
-        * x        (x) <U1 'a' 'b' 'c'
+          * x        (x) <U1 'a' 'b' 'c'
         >>> array
         <xarray.DataArray (x: 3)>
         array([7, 2, 3])
         Coordinates:
-        * x        (x) <U1 'a' 'b' 'c'
+          * x        (x) <U1 'a' 'b' 'c'
 
         Changing the data using the ``data`` argument maintains the
         structure of the original object, but with the new data. Original
@@ -908,14 +908,14 @@ def copy(self, deep: bool = True, data: Any = None) -> "DataArray":
 
         >>> array.copy(data=[0.1, 0.2, 0.3])
         <xarray.DataArray (x: 3)>
-        array([ 0.1,  0.2,  0.3])
+        array([0.1, 0.2, 0.3])
         Coordinates:
-        * x        (x) <U1 'a' 'b' 'c'
+          * x        (x) <U1 'a' 'b' 'c'
         >>> array
         <xarray.DataArray (x: 3)>
-        array([1, 2, 3])
+        array([7, 2, 3])
         Coordinates:
-        * x        (x) <U1 'a' 'b' 'c'
+          * x        (x) <U1 'a' 'b' 'c'
 
         See Also
         --------
@@ -1231,26 +1231,36 @@ def broadcast_like(
         Examples
         --------
 
+        >>> arr1 = xr.DataArray(
+        ...     np.random.randn(2, 3),
+        ...     dims=("x", "y"),
+        ...     coords={"x": ["a", "b"], "y": ["a", "b", "c"]},
+        ... )
+        >>> arr2 = xr.DataArray(
+        ...     np.random.randn(3, 2),
+        ...     dims=("x", "y"),
+        ...     coords={"x": ["a", "b", "c"], "y": ["a", "b"]},
+        ... )
         >>> arr1
         <xarray.DataArray (x: 2, y: 3)>
-        array([[0.840235, 0.215216, 0.77917 ],
-               [0.726351, 0.543824, 0.875115]])
+        array([[ 1.76405235,  0.40015721,  0.97873798],
+               [ 2.2408932 ,  1.86755799, -0.97727788]])
         Coordinates:
           * x        (x) <U1 'a' 'b'
           * y        (y) <U1 'a' 'b' 'c'
         >>> arr2
         <xarray.DataArray (x: 3, y: 2)>
-        array([[0.612611, 0.125753],
-               [0.853181, 0.948818],
-               [0.180885, 0.33363 ]])
+        array([[ 0.95008842, -0.15135721],
+               [-0.10321885,  0.4105985 ],
+               [ 0.14404357,  1.45427351]])
         Coordinates:
           * x        (x) <U1 'a' 'b' 'c'
           * y        (y) <U1 'a' 'b'
         >>> arr1.broadcast_like(arr2)
         <xarray.DataArray (x: 3, y: 3)>
-        array([[0.840235, 0.215216, 0.77917 ],
-               [0.726351, 0.543824, 0.875115],
-               [     nan,      nan,      nan]])
+        array([[ 1.76405235,  0.40015721,  0.97873798],
+               [ 2.2408932 ,  1.86755799, -0.97727788],
+               [        nan,         nan,         nan]])
         Coordinates:
           * x        (x) object 'a' 'b' 'c'
           * y        (y) object 'a' 'b' 'c'
@@ -1453,7 +1463,7 @@ def interp(
         >>> da = xr.DataArray([1, 3], [("x", np.arange(2))])
         >>> da.interp(x=0.5)
         <xarray.DataArray ()>
-        array(2.0)
+        array(2.)
         Coordinates:
             x        float64 0.5
         """
@@ -1584,7 +1594,9 @@ def swap_dims(self, dims_dict: Mapping[Hashable, Hashable]) -> "DataArray":
         --------
 
         >>> arr = xr.DataArray(
-        ...     data=[0, 1], dims="x", coords={"x": ["a", "b"], "y": ("x", [0, 1])},
+        ...     data=[0, 1],
+        ...     dims="x",
+        ...     coords={"x": ["a", "b"], "y": ("x", [0, 1])},
         ... )
         >>> arr
         <xarray.DataArray (x: 2)>
@@ -1840,12 +1852,16 @@ def stack(
         array([[0, 1, 2],
                [3, 4, 5]])
         Coordinates:
-          * x        (x) |S1 'a' 'b'
+          * x        (x) <U1 'a' 'b'
           * y        (y) int64 0 1 2
         >>> stacked = arr.stack(z=("x", "y"))
         >>> stacked.indexes["z"]
-        MultiIndex(levels=[['a', 'b'], [0, 1, 2]],
-                   codes=[[0, 0, 0, 1, 1, 1], [0, 1, 2, 0, 1, 2]],
+        MultiIndex([('a', 0),
+                    ('a', 1),
+                    ('a', 2),
+                    ('b', 0),
+                    ('b', 1),
+                    ('b', 2)],
                    names=['x', 'y'])
 
         See Also
@@ -1897,12 +1913,16 @@ def unstack(
         array([[0, 1, 2],
                [3, 4, 5]])
         Coordinates:
-          * x        (x) |S1 'a' 'b'
+          * x        (x) <U1 'a' 'b'
           * y        (y) int64 0 1 2
         >>> stacked = arr.stack(z=("x", "y"))
         >>> stacked.indexes["z"]
-        MultiIndex(levels=[['a', 'b'], [0, 1, 2]],
-                   codes=[[0, 0, 0, 1, 1, 1], [0, 1, 2, 0, 1, 2]],
+        MultiIndex([('a', 0),
+                    ('a', 1),
+                    ('a', 2),
+                    ('b', 0),
+                    ('b', 1),
+                    ('b', 2)],
                    names=['x', 'y'])
         >>> roundtripped = stacked.unstack()
         >>> arr.identical(roundtripped)
@@ -1953,11 +1973,13 @@ def to_unstacked_dataset(self, dim, level=0):
         Data variables:
             a        (x, y) int64 0 1 2 3 4 5
             b        (x) int64 0 3
-        >>> stacked = data.to_stacked_array("z", ["y"])
+        >>> stacked = data.to_stacked_array("z", ["x"])
         >>> stacked.indexes["z"]
-        MultiIndex(levels=[['a', 'b'], [0, 1, 2]],
-                labels=[[0, 0, 0, 1], [0, 1, 2, -1]],
-                names=['variable', 'y'])
+        MultiIndex([('a', 0.0),
+                    ('a', 1.0),
+                    ('a', 2.0),
+                    ('b', nan)],
+                   names=['variable', 'y'])
         >>> roundtripped = stacked.to_unstacked_dataset(dim="z")
         >>> data.identical(roundtripped)
         True
@@ -2842,12 +2864,12 @@ def diff(self, dim: Hashable, n: int = 1, label: Hashable = "upper") -> "DataArr
         <xarray.DataArray (x: 3)>
         array([0, 1, 0])
         Coordinates:
-        * x        (x) int64 2 3 4
+          * x        (x) int64 2 3 4
         >>> arr.diff("x", 2)
         <xarray.DataArray (x: 2)>
         array([ 1, -1])
         Coordinates:
-        * x        (x) int64 3 4
+          * x        (x) int64 3 4
 
         See Also
         --------
@@ -2896,9 +2918,8 @@ def shift(
         >>> arr = xr.DataArray([5, 6, 7], dims="x")
         >>> arr.shift(x=1)
         <xarray.DataArray (x: 3)>
-        array([ nan,   5.,   6.])
-        Coordinates:
-          * x        (x) int64 0 1 2
+        array([nan,  5.,  6.])
+        Dimensions without coordinates: x
         """
         variable = self.variable.shift(
             shifts=shifts, fill_value=fill_value, **shifts_kwargs
@@ -2948,8 +2969,7 @@ def roll(
         >>> arr.roll(x=1)
         <xarray.DataArray (x: 3)>
         array([7, 5, 6])
-        Coordinates:
-          * x        (x) int64 2 0 1
+        Dimensions without coordinates: x
         """
         ds = self._to_temp_dataset().roll(
             shifts=shifts, roll_coords=roll_coords, **shifts_kwargs
@@ -2998,7 +3018,7 @@ def dot(
         >>> dm = xr.DataArray(dm_vals, dims=["z"])
 
         >>> dm.dims
-        ('z')
+        ('z',)
 
         >>> da.dims
         ('x', 'y', 'z')
@@ -3062,15 +3082,15 @@ def sortby(
         ... )
         >>> da
         <xarray.DataArray (time: 5)>
-        array([ 0.965471,  0.615637,  0.26532 ,  0.270962,  0.552878])
+        array([0.5488135 , 0.71518937, 0.60276338, 0.54488318, 0.4236548 ])
         Coordinates:
-          * time     (time) datetime64[ns] 2000-01-01 2000-01-02 2000-01-03 ...
+          * time     (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2000-01-05
 
         >>> da.sortby(da)
         <xarray.DataArray (time: 5)>
-        array([ 0.26532 ,  0.270962,  0.552878,  0.615637,  0.965471])
+        array([0.4236548 , 0.54488318, 0.5488135 , 0.60276338, 0.71518937])
         Coordinates:
-          * time     (time) datetime64[ns] 2000-01-03 2000-01-04 2000-01-05 ...
+          * time     (time) datetime64[ns] 2000-01-05 2000-01-04 ... 2000-01-02
         """
         ds = self._to_temp_dataset().sortby(variables, ascending=ascending)
         return self._from_temp_dataset(ds)
@@ -3203,7 +3223,7 @@ def rank(
         >>> arr = xr.DataArray([5, 6, 7], dims="x")
         >>> arr.rank("x")
         <xarray.DataArray (x: 3)>
-        array([ 1.,   2.,   3.])
+        array([1., 2., 3.])
         Dimensions without coordinates: x
         """
 
@@ -3258,10 +3278,10 @@ def differentiate(
         >>>
         >>> da.differentiate("x")
         <xarray.DataArray (x: 4, y: 3)>
-        array([[30.      , 30.      , 30.      ],
-               [27.545455, 27.545455, 27.545455],
-               [27.545455, 27.545455, 27.545455],
-               [30.      , 30.      , 30.      ]])
+        array([[30.        , 30.        , 30.        ],
+               [27.54545455, 27.54545455, 27.54545455],
+               [27.54545455, 27.54545455, 27.54545455],
+               [30.        , 30.        , 30.        ]])
         Coordinates:
           * x        (x) float64 0.0 0.1 1.1 1.2
         Dimensions without coordinates: y
@@ -3428,7 +3448,7 @@ def map_blocks(
         to the function being applied in ``xr.map_blocks()``:
 
         >>> array.map_blocks(
-        ...     calculate_anomaly, kwargs={"groupby_type": "time.year"}, template=array,
+        ...     calculate_anomaly, kwargs={"groupby_type": "time.year"}, template=array
         ... )  # doctest: +ELLIPSIS
         <xarray.DataArray (time: 24)>
         dask.array<calculate_anomaly-...-<this, shape=(24,), dtype=float64, chunksize=(24,), chunktype=numpy.ndarray>
@@ -3924,6 +3944,7 @@ def argmin(
         {'x': <xarray.DataArray ()>
         array(2)}
         >>> array.isel(array.argmin(...))
+        <xarray.DataArray ()>
         array(-1)
 
         >>> array = xr.DataArray(
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 92de628f5ad..825d2044a12 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -1018,10 +1018,10 @@ def copy(self, deep: bool = False, data: Mapping = None) -> "Dataset":
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
         Coordinates:
-        * x        (x) <U3 'one' 'two'
+          * x        (x) <U3 'one' 'two'
         Dimensions without coordinates: dim_0, dim_1
         Data variables:
-            foo      (dim_0, dim_1) float64 -0.8079 0.3897 -1.862 -0.6091 -1.051 -0.3003
+            foo      (dim_0, dim_1) float64 1.764 0.4002 0.9787 2.241 1.868 -0.9773
             bar      (x) int64 -1 2
 
         >>> ds_0 = ds.copy(deep=False)
@@ -1030,20 +1030,20 @@ def copy(self, deep: bool = False, data: Mapping = None) -> "Dataset":
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
         Coordinates:
-        * x        (x) <U3 'one' 'two'
+          * x        (x) <U3 'one' 'two'
         Dimensions without coordinates: dim_0, dim_1
         Data variables:
-            foo      (dim_0, dim_1) float64 7.0 0.3897 -1.862 -0.6091 -1.051 -0.3003
+            foo      (dim_0, dim_1) float64 7.0 0.4002 0.9787 2.241 1.868 -0.9773
             bar      (x) int64 -1 2
 
         >>> ds
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
         Coordinates:
-        * x        (x) <U3 'one' 'two'
+          * x        (x) <U3 'one' 'two'
         Dimensions without coordinates: dim_0, dim_1
         Data variables:
-            foo      (dim_0, dim_1) float64 7.0 0.3897 -1.862 -0.6091 -1.051 -0.3003
+            foo      (dim_0, dim_1) float64 7.0 0.4002 0.9787 2.241 1.868 -0.9773
             bar      (x) int64 -1 2
 
         Changing the data using the ``data`` argument maintains the
@@ -1054,7 +1054,7 @@ def copy(self, deep: bool = False, data: Mapping = None) -> "Dataset":
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
         Coordinates:
-        * x        (x) <U3 'one' 'two'
+          * x        (x) <U3 'one' 'two'
         Dimensions without coordinates: dim_0, dim_1
         Data variables:
             foo      (dim_0, dim_1) int64 0 1 2 3 4 5
@@ -1064,10 +1064,10 @@ def copy(self, deep: bool = False, data: Mapping = None) -> "Dataset":
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
         Coordinates:
-        * x        (x) <U3 'one' 'two'
+          * x        (x) <U3 'one' 'two'
         Dimensions without coordinates: dim_0, dim_1
         Data variables:
-            foo      (dim_0, dim_1) float64 7.0 0.3897 -1.862 -0.6091 -1.051 -0.3003
+            foo      (dim_0, dim_1) float64 7.0 0.4002 0.9787 2.241 1.868 -0.9773
             bar      (x) int64 -1 2
 
         See Also
@@ -2402,10 +2402,10 @@ def reindex(
         <xarray.Dataset>
         Dimensions:      (station: 4)
         Coordinates:
-        * station      (station) <U7 'boston' 'nyc' 'seattle' 'denver'
+          * station      (station) <U7 'boston' 'nyc' 'seattle' 'denver'
         Data variables:
-            temperature  (station) float64 18.84 14.59 19.22 17.16
-            pressure     (station) float64 324.1 194.3 122.8 244.3
+            temperature  (station) float64 10.98 14.3 12.06 10.9
+            pressure     (station) float64 211.8 322.9 218.8 445.9
         >>> x.indexes
         station: Index(['boston', 'nyc', 'seattle', 'denver'], dtype='object', name='station')
 
@@ -2417,10 +2417,10 @@ def reindex(
         <xarray.Dataset>
         Dimensions:      (station: 4)
         Coordinates:
-        * station      (station) object 'boston' 'austin' 'seattle' 'lincoln'
+          * station      (station) object 'boston' 'austin' 'seattle' 'lincoln'
         Data variables:
-            temperature  (station) float64 18.84 nan 19.22 nan
-            pressure     (station) float64 324.1 nan 122.8 nan
+            temperature  (station) float64 10.98 nan 12.06 nan
+            pressure     (station) float64 211.8 nan 218.8 nan
 
         We can fill in the missing values by passing a value to the keyword `fill_value`.
 
@@ -2428,10 +2428,10 @@ def reindex(
         <xarray.Dataset>
         Dimensions:      (station: 4)
         Coordinates:
-        * station      (station) object 'boston' 'austin' 'seattle' 'lincoln'
+          * station      (station) object 'boston' 'austin' 'seattle' 'lincoln'
         Data variables:
-            temperature  (station) float64 18.84 0.0 19.22 0.0
-            pressure     (station) float64 324.1 0.0 122.8 0.0
+            temperature  (station) float64 10.98 0.0 12.06 0.0
+            pressure     (station) float64 211.8 0.0 218.8 0.0
 
         We can also use different fill values for each variable.
 
@@ -2441,10 +2441,10 @@ def reindex(
         <xarray.Dataset>
         Dimensions:      (station: 4)
         Coordinates:
-        * station      (station) object 'boston' 'austin' 'seattle' 'lincoln'
+          * station      (station) object 'boston' 'austin' 'seattle' 'lincoln'
         Data variables:
-            temperature  (station) float64 18.84 0.0 19.22 0.0
-            pressure     (station) float64 324.1 100.0 122.8 100.0
+            temperature  (station) float64 10.98 0.0 12.06 0.0
+            pressure     (station) float64 211.8 100.0 218.8 100.0
 
         Because the index is not monotonically increasing or decreasing, we cannot use arguments
         to the keyword method to fill the `NaN` values.
@@ -2472,10 +2472,10 @@ def reindex(
         <xarray.Dataset>
         Dimensions:      (time: 6)
         Coordinates:
-        * time         (time) datetime64[ns] 2019-01-01 2019-01-02 ... 2019-01-06
+          * time         (time) datetime64[ns] 2019-01-01 2019-01-02 ... 2019-01-06
         Data variables:
             temperature  (time) float64 15.57 12.77 nan 0.3081 16.59 15.12
-            pressure     (time) float64 103.4 122.7 452.0 444.0 399.2 486.0
+            pressure     (time) float64 481.8 191.7 395.9 264.4 284.0 462.8
 
         Suppose we decide to expand the dataset to cover a wider date range.
 
@@ -2484,10 +2484,10 @@ def reindex(
         <xarray.Dataset>
         Dimensions:      (time: 10)
         Coordinates:
-        * time         (time) datetime64[ns] 2018-12-29 2018-12-30 ... 2019-01-07
+          * time         (time) datetime64[ns] 2018-12-29 2018-12-30 ... 2019-01-07
         Data variables:
             temperature  (time) float64 nan nan nan 15.57 ... 0.3081 16.59 15.12 nan
-            pressure     (time) float64 nan nan nan 103.4 ... 444.0 399.2 486.0 nan
+            pressure     (time) float64 nan nan nan 481.8 ... 264.4 284.0 462.8 nan
 
         The index entries that did not have a value in the original data frame (for example, `2018-12-29`)
         are by default filled with NaN. If desired, we can fill in the missing values using one of several options.
@@ -2500,10 +2500,10 @@ def reindex(
         <xarray.Dataset>
         Dimensions:      (time: 10)
         Coordinates:
-        * time         (time) datetime64[ns] 2018-12-29 2018-12-30 ... 2019-01-07
+          * time         (time) datetime64[ns] 2018-12-29 2018-12-30 ... 2019-01-07
         Data variables:
             temperature  (time) float64 15.57 15.57 15.57 15.57 ... 16.59 15.12 nan
-            pressure     (time) float64 103.4 103.4 103.4 103.4 ... 399.2 486.0 nan
+            pressure     (time) float64 481.8 481.8 481.8 481.8 ... 284.0 462.8 nan
 
         Please note that the `NaN` value present in the original dataset (at index value `2019-01-03`)
         will not be filled by any of the value propagation schemes.
@@ -2512,18 +2512,18 @@ def reindex(
         <xarray.Dataset>
         Dimensions:      (time: 1)
         Coordinates:
-        * time         (time) datetime64[ns] 2019-01-03
+          * time         (time) datetime64[ns] 2019-01-03
         Data variables:
             temperature  (time) float64 nan
-            pressure     (time) float64 452.0
+            pressure     (time) float64 395.9
         >>> x3.where(x3.temperature.isnull(), drop=True)
         <xarray.Dataset>
         Dimensions:      (time: 2)
         Coordinates:
-        * time         (time) datetime64[ns] 2019-01-03 2019-01-07
+          * time         (time) datetime64[ns] 2019-01-03 2019-01-07
         Data variables:
             temperature  (time) float64 nan nan
-            pressure     (time) float64 452.0 nan
+            pressure     (time) float64 395.9 nan
 
         This is because filling while reindexing does not look at dataset values, but only compares
         the original and desired indexes. If you do want to fill in the `NaN` values present in the
@@ -3431,20 +3431,20 @@ def to_stacked_array(
         <xarray.Dataset>
         Dimensions:  (x: 2, y: 3)
         Coordinates:
-        * y        (y) <U1 'u' 'v' 'w'
+          * y        (y) <U1 'u' 'v' 'w'
         Dimensions without coordinates: x
         Data variables:
             a        (x, y) int64 0 1 2 3 4 5
             b        (x) int64 6 7
 
         >>> data.to_stacked_array("z", sample_dims=["x"])
-        <xarray.DataArray (x: 2, z: 4)>
+        <xarray.DataArray 'a' (x: 2, z: 4)>
         array([[0, 1, 2, 6],
-            [3, 4, 5, 7]])
+               [3, 4, 5, 7]])
         Coordinates:
-        * z         (z) MultiIndex
-        - variable  (z) object 'a' 'a' 'a' 'b'
-        - y         (z) object 'u' 'v' 'w' nan
+          * z         (z) MultiIndex
+          - variable  (z) object 'a' 'a' 'a' 'b'
+          - y         (z) object 'u' 'v' 'w' nan
         Dimensions without coordinates: x
 
         """
@@ -3822,7 +3822,7 @@ def drop_sel(self, labels=None, *, errors="raise", **labels_kwargs):
           * y        (y) <U1 'b'
         Dimensions without coordinates: x
         Data variables:
-            A        (x, y) float64 -0.3454 0.1734
+            A        (x, y) float64 0.4002 1.868
         >>> ds.drop_sel(y="b")
         <xarray.Dataset>
         Dimensions:  (x: 2, y: 2)
@@ -3830,7 +3830,7 @@ def drop_sel(self, labels=None, *, errors="raise", **labels_kwargs):
           * y        (y) <U1 'a' 'c'
         Dimensions without coordinates: x
         Data variables:
-            A        (x, y) float64 -0.3944 -1.418 1.423 -1.041
+            A        (x, y) float64 1.764 0.9787 2.241 -0.9773
         """
         if errors not in ["raise", "ignore"]:
             raise ValueError('errors must be either "raise" or "ignore"')
@@ -4036,7 +4036,7 @@ def fillna(self, value: Any) -> "Dataset":
         <xarray.Dataset>
         Dimensions:  (x: 4)
         Coordinates:
-        * x        (x) int64 0 1 2 3
+          * x        (x) int64 0 1 2 3
         Data variables:
             A        (x) float64 nan 2.0 nan 0.0
             B        (x) float64 3.0 4.0 nan 1.0
@@ -4049,7 +4049,7 @@ def fillna(self, value: Any) -> "Dataset":
         <xarray.Dataset>
         Dimensions:  (x: 4)
         Coordinates:
-        * x        (x) int64 0 1 2 3
+          * x        (x) int64 0 1 2 3
         Data variables:
             A        (x) float64 0.0 2.0 0.0 0.0
             B        (x) float64 3.0 4.0 0.0 1.0
@@ -4063,7 +4063,7 @@ def fillna(self, value: Any) -> "Dataset":
         <xarray.Dataset>
         Dimensions:  (x: 4)
         Coordinates:
-        * x        (x) int64 0 1 2 3
+          * x        (x) int64 0 1 2 3
         Data variables:
             A        (x) float64 0.0 2.0 0.0 0.0
             B        (x) float64 3.0 4.0 1.0 1.0
@@ -4371,14 +4371,14 @@ def map(
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
         Dimensions without coordinates: dim_0, dim_1, x
         Data variables:
-            foo      (dim_0, dim_1) float64 -0.3751 -1.951 -1.945 0.2948 0.711 -0.3948
+            foo      (dim_0, dim_1) float64 1.764 0.4002 0.9787 2.241 1.868 -0.9773
             bar      (x) int64 -1 2
         >>> ds.map(np.fabs)
         <xarray.Dataset>
         Dimensions:  (dim_0: 2, dim_1: 3, x: 2)
         Dimensions without coordinates: dim_0, dim_1, x
         Data variables:
-            foo      (dim_0, dim_1) float64 0.3751 1.951 1.945 0.2948 0.711 0.3948
+            foo      (dim_0, dim_1) float64 1.764 0.4002 0.9787 2.241 1.868 0.9773
             bar      (x) float64 1.0 2.0
         """
         variables = {
@@ -4462,11 +4462,11 @@ def assign(
         <xarray.Dataset>
         Dimensions:        (lat: 2, lon: 2)
         Coordinates:
-        * lat            (lat) int64 10 20
-        * lon            (lon) int64 150 160
+          * lat            (lat) int64 10 20
+          * lon            (lon) int64 150 160
         Data variables:
-            temperature_c  (lat, lon) float64 18.04 12.51 17.64 9.313
-            precipitation  (lat, lon) float64 0.4751 0.6827 0.3697 0.03524
+            temperature_c  (lat, lon) float64 10.98 14.3 12.06 10.9
+            precipitation  (lat, lon) float64 0.4237 0.6459 0.4376 0.8918
 
         Where the value is a callable, evaluated on dataset:
 
@@ -4474,12 +4474,12 @@ def assign(
         <xarray.Dataset>
         Dimensions:        (lat: 2, lon: 2)
         Coordinates:
-        * lat            (lat) int64 10 20
-        * lon            (lon) int64 150 160
+          * lat            (lat) int64 10 20
+          * lon            (lon) int64 150 160
         Data variables:
-            temperature_c  (lat, lon) float64 18.04 12.51 17.64 9.313
-            precipitation  (lat, lon) float64 0.4751 0.6827 0.3697 0.03524
-            temperature_f  (lat, lon) float64 64.47 54.51 63.75 48.76
+            temperature_c  (lat, lon) float64 10.98 14.3 12.06 10.9
+            precipitation  (lat, lon) float64 0.4237 0.6459 0.4376 0.8918
+            temperature_f  (lat, lon) float64 51.76 57.75 53.7 51.62
 
         Alternatively, the same behavior can be achieved by directly referencing an existing dataarray:
 
@@ -4487,12 +4487,12 @@ def assign(
         <xarray.Dataset>
         Dimensions:        (lat: 2, lon: 2)
         Coordinates:
-        * lat            (lat) int64 10 20
-        * lon            (lon) int64 150 160
+          * lat            (lat) int64 10 20
+          * lon            (lon) int64 150 160
         Data variables:
-            temperature_c  (lat, lon) float64 18.04 12.51 17.64 9.313
-            precipitation  (lat, lon) float64 0.4751 0.6827 0.3697 0.03524
-            temperature_f  (lat, lon) float64 64.47 54.51 63.75 48.76
+            temperature_c  (lat, lon) float64 10.98 14.3 12.06 10.9
+            precipitation  (lat, lon) float64 0.4237 0.6459 0.4376 0.8918
+            temperature_f  (lat, lon) float64 51.76 57.75 53.7 51.62
 
         """
         variables = either_dict_or_kwargs(variables, variables_kwargs, "assign")
@@ -5056,17 +5056,15 @@ def diff(self, dim, n=1, label="upper"):
         >>> ds.diff("x")
         <xarray.Dataset>
         Dimensions:  (x: 3)
-        Coordinates:
-          * x        (x) int64 1 2 3
+        Dimensions without coordinates: x
         Data variables:
             foo      (x) int64 0 1 0
         >>> ds.diff("x", 2)
         <xarray.Dataset>
         Dimensions:  (x: 2)
-        Coordinates:
-        * x        (x) int64 2 3
+        Dimensions without coordinates: x
         Data variables:
-        foo      (x) int64 1 -1
+            foo      (x) int64 1 -1
 
         See Also
         --------
@@ -5149,8 +5147,7 @@ def shift(self, shifts=None, fill_value=dtypes.NA, **shifts_kwargs):
         >>> ds.shift(x=2)
         <xarray.Dataset>
         Dimensions:  (x: 5)
-        Coordinates:
-          * x        (x) int64 0 1 2 3 4
+        Dimensions without coordinates: x
         Data variables:
             foo      (x) object nan nan 'a' 'b' 'c'
         """
@@ -5214,10 +5211,9 @@ def roll(self, shifts=None, roll_coords=None, **shifts_kwargs):
         >>> ds.roll(x=2)
         <xarray.Dataset>
         Dimensions:  (x: 5)
-        Coordinates:
-          * x        (x) int64 3 4 0 1 2
+        Dimensions without coordinates: x
         Data variables:
-            foo      (x) object 'd' 'e' 'a' 'b' 'c'
+            foo      (x) <U1 'd' 'e' 'a' 'b' 'c'
         """
         shifts = either_dict_or_kwargs(shifts, shifts_kwargs, "roll")
         invalid = [k for k in shifts if k not in self.dims]
@@ -5733,29 +5729,27 @@ def filter_by_attrs(self, **kwargs):
         <xarray.Dataset>
         Dimensions:         (time: 3, x: 2, y: 2)
         Coordinates:
-          * x               (x) int64 0 1
-          * time            (time) datetime64[ns] 2014-09-06 2014-09-07 2014-09-08
-            lat             (x, y) float64 42.25 42.21 42.63 42.59
-          * y               (y) int64 0 1
             reference_time  datetime64[ns] 2014-09-05
+            lat             (x, y) float64 42.25 42.21 42.63 42.59
+          * time            (time) datetime64[ns] 2014-09-06 2014-09-07 2014-09-08
             lon             (x, y) float64 -99.83 -99.32 -99.79 -99.23
+        Dimensions without coordinates: x, y
         Data variables:
-            precipitation   (x, y, time) float64 4.178 2.307 6.041 6.046 0.06648 ...
+            precipitation   (x, y, time) float64 5.68 9.256 0.7104 ... 7.992 4.615 7.805
         >>> # Get all variables that have a standard_name attribute.
         >>> standard_name = lambda v: v is not None
         >>> ds.filter_by_attrs(standard_name=standard_name)
         <xarray.Dataset>
         Dimensions:         (time: 3, x: 2, y: 2)
         Coordinates:
-            lon             (x, y) float64 -99.83 -99.32 -99.79 -99.23
+            reference_time  datetime64[ns] 2014-09-05
             lat             (x, y) float64 42.25 42.21 42.63 42.59
-          * x               (x) int64 0 1
-          * y               (y) int64 0 1
           * time            (time) datetime64[ns] 2014-09-06 2014-09-07 2014-09-08
-            reference_time  datetime64[ns] 2014-09-05
+            lon             (x, y) float64 -99.83 -99.32 -99.79 -99.23
+        Dimensions without coordinates: x, y
         Data variables:
-            temperature     (x, y, time) float64 25.86 20.82 6.954 23.13 10.25 11.68 ...
-            precipitation   (x, y, time) float64 5.702 0.9422 2.075 1.178 3.284 ...
+            temperature     (x, y, time) float64 29.11 18.2 22.83 ... 18.28 16.15 26.63
+            precipitation   (x, y, time) float64 5.68 9.256 0.7104 ... 7.992 4.615 7.805
 
         """
         selection = []
@@ -6370,7 +6364,7 @@ def idxmin(
           * y        (y) int64 -1 0 1
         Data variables:
             int      <U1 'e'
-            float    (y) <U1 'e' 'a' 'c'
+            float    (y) object 'e' 'a' 'c'
         """
         return self.map(
             methodcaller(
diff --git a/xarray/core/extensions.py b/xarray/core/extensions.py
index e81070d18fd..ee4c3ebc9e6 100644
--- a/xarray/core/extensions.py
+++ b/xarray/core/extensions.py
@@ -88,35 +88,32 @@ def register_dataset_accessor(name):
     Examples
     --------
 
-    In your library code::
-
-        import xarray as xr
-
-        @xr.register_dataset_accessor('geo')
-        class GeoAccessor:
-            def __init__(self, xarray_obj):
-                self._obj = xarray_obj
-
-            @property
-            def center(self):
-                # return the geographic center point of this dataset
-                lon = self._obj.latitude
-                lat = self._obj.longitude
-                return (float(lon.mean()), float(lat.mean()))
-
-            def plot(self):
-                # plot this array's data on a map, e.g., using Cartopy
-                pass
+    In your library code:
+
+    >>> @xr.register_dataset_accessor("geo")
+    ... class GeoAccessor:
+    ...     def __init__(self, xarray_obj):
+    ...         self._obj = xarray_obj
+    ...
+    ...     @property
+    ...     def center(self):
+    ...         # return the geographic center point of this dataset
+    ...         lon = self._obj.latitude
+    ...         lat = self._obj.longitude
+    ...         return (float(lon.mean()), float(lat.mean()))
+    ...
+    ...     def plot(self):
+    ...         # plot this array's data on a map, e.g., using Cartopy
+    ...         pass
 
     Back in an interactive IPython session:
 
-        >>> ds = xarray.Dataset(
-        ...     {"longitude": np.linspace(0, 10), "latitude": np.linspace(0, 20)}
-        ... )
-        >>> ds.geo.center
-        (5.0, 10.0)
-        >>> ds.geo.plot()
-        # plots data on a map
+    >>> ds = xr.Dataset(
+    ...     {"longitude": np.linspace(0, 10), "latitude": np.linspace(0, 20)}
+    ... )
+    >>> ds.geo.center
+    (10.0, 5.0)
+    >>> ds.geo.plot()  # plots data on a map
 
     See also
     --------
diff --git a/xarray/core/groupby.py b/xarray/core/groupby.py
index a5d96bc66cc..5f328d7a03a 100644
--- a/xarray/core/groupby.py
+++ b/xarray/core/groupby.py
@@ -599,7 +599,7 @@ def quantile(
         >>> da = xr.DataArray(
         ...     [[1.3, 8.4, 0.7, 6.9], [0.7, 4.2, 9.4, 1.5], [6.5, 7.3, 2.6, 1.9]],
         ...     coords={"x": [0, 0, 1], "y": [1, 1, 2, 2]},
-        ...     dims=("y", "y"),
+        ...     dims=("x", "y"),
         ... )
         >>> ds = xr.Dataset({"a": da})
         >>> da.groupby("x").quantile(0)
@@ -624,13 +624,14 @@ def quantile(
                 [4.2 , 6.3 , 8.4 ],
                 [0.7 , 5.05, 9.4 ],
                 [1.5 , 4.2 , 6.9 ]],
+        <BLANKLINE>
                [[6.5 , 6.5 , 6.5 ],
                 [7.3 , 7.3 , 7.3 ],
                 [2.6 , 2.6 , 2.6 ],
                 [1.9 , 1.9 , 1.9 ]]])
         Coordinates:
-          * y         (y) int64 1 1 2 2
           * quantile  (quantile) float64 0.0 0.5 1.0
+          * y         (y) int64 1 1 2 2
           * x         (x) int64 0 1
         >>> ds.groupby("y").quantile([0, 0.5, 1], dim=...)
         <xarray.Dataset>
diff --git a/xarray/core/indexing.py b/xarray/core/indexing.py
index 919a9cf5293..66c62653139 100644
--- a/xarray/core/indexing.py
+++ b/xarray/core/indexing.py
@@ -899,10 +899,14 @@ def _decompose_vectorized_indexer(
     Even if the backend array only supports outer indexing, it is more
     efficient to load a subslice of the array than loading the entire array,
 
-    >>> backend_indexer = OuterIndexer([0, 1, 3], [2, 3])
-    >>> array = array[backend_indexer]  # load subslice of the array
-    >>> np_indexer = VectorizedIndexer([0, 2, 1], [0, 1, 0])
-    >>> array[np_indexer]  # vectorized indexing for on-memory np.ndarray.
+    >>> array = np.arange(36).reshape(6, 6)
+    >>> backend_indexer = OuterIndexer((np.array([0, 1, 3]), np.array([2, 3])))
+    >>> # load subslice of the array
+    ... array = NumpyIndexingAdapter(array)[backend_indexer]
+    >>> np_indexer = VectorizedIndexer((np.array([0, 2, 1]), np.array([0, 1, 0])))
+    >>> # vectorized indexing for on-memory np.ndarray.
+    ... NumpyIndexingAdapter(array)[np_indexer]
+    array([ 2, 21,  8])
     """
     assert isinstance(indexer, VectorizedIndexer)
 
@@ -977,10 +981,16 @@ def _decompose_outer_indexer(
     Even if the backend array only supports basic indexing, it is more
     efficient to load a subslice of the array than loading the entire array,
 
-    >>> backend_indexer = BasicIndexer(slice(0, 3), slice(2, 3))
-    >>> array = array[backend_indexer]  # load subslice of the array
-    >>> np_indexer = OuterIndexer([0, 2, 1], [0, 1, 0])
-    >>> array[np_indexer]  # outer indexing for on-memory np.ndarray.
+    >>> array = np.arange(36).reshape(6, 6)
+    >>> backend_indexer = BasicIndexer((slice(0, 3), slice(2, 4)))
+    >>> # load subslice of the array
+    ... array = NumpyIndexingAdapter(array)[backend_indexer]
+    >>> np_indexer = OuterIndexer((np.array([0, 2, 1]), np.array([0, 1, 0])))
+    >>> # outer indexing for on-memory np.ndarray.
+    ... NumpyIndexingAdapter(array)[np_indexer]
+    array([[ 2,  3,  2],
+           [14, 15, 14],
+           [ 8,  9,  8]])
     """
     if indexing_support == IndexingSupport.VECTORIZED:
         return indexer, BasicIndexer(())
diff --git a/xarray/core/merge.py b/xarray/core/merge.py
index 231e1f7db10..ca4e29b600d 100644
--- a/xarray/core/merge.py
+++ b/xarray/core/merge.py
@@ -711,32 +711,32 @@ def merge(
     array([[1., 2.],
            [3., 5.]])
     Coordinates:
-    * lat      (lat) float64 35.0 40.0
-    * lon      (lon) float64 100.0 120.0
+      * lat      (lat) float64 35.0 40.0
+      * lon      (lon) float64 100.0 120.0
 
     >>> y
     <xarray.DataArray 'var2' (lat: 2, lon: 2)>
     array([[5., 6.],
            [7., 8.]])
     Coordinates:
-    * lat      (lat) float64 35.0 42.0
-    * lon      (lon) float64 100.0 150.0
+      * lat      (lat) float64 35.0 42.0
+      * lon      (lon) float64 100.0 150.0
 
     >>> z
     <xarray.DataArray 'var3' (time: 2, lon: 2)>
     array([[0., 3.],
            [4., 9.]])
     Coordinates:
-    * time     (time) float64 30.0 60.0
-    * lon      (lon) float64 100.0 150.0
+      * time     (time) float64 30.0 60.0
+      * lon      (lon) float64 100.0 150.0
 
     >>> xr.merge([x, y, z])
     <xarray.Dataset>
     Dimensions:  (lat: 3, lon: 3, time: 2)
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0 150.0
-    * time     (time) float64 30.0 60.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0 150.0
+      * time     (time) float64 30.0 60.0
     Data variables:
         var1     (lat, lon) float64 1.0 2.0 nan 3.0 5.0 nan nan nan nan
         var2     (lat, lon) float64 5.0 nan 6.0 nan nan nan 7.0 nan 8.0
@@ -746,9 +746,9 @@ def merge(
     <xarray.Dataset>
     Dimensions:  (lat: 3, lon: 3, time: 2)
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0 150.0
-    * time     (time) float64 30.0 60.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0 150.0
+      * time     (time) float64 30.0 60.0
     Data variables:
         var1     (lat, lon) float64 1.0 2.0 nan 3.0 5.0 nan nan nan nan
         var2     (lat, lon) float64 5.0 nan 6.0 nan nan nan 7.0 nan 8.0
@@ -758,9 +758,9 @@ def merge(
     <xarray.Dataset>
     Dimensions:  (lat: 3, lon: 3, time: 2)
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0 150.0
-    * time     (time) float64 30.0 60.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0 150.0
+      * time     (time) float64 30.0 60.0
     Data variables:
         var1     (lat, lon) float64 1.0 2.0 nan 3.0 5.0 nan nan nan nan
         var2     (lat, lon) float64 5.0 nan 6.0 nan nan nan 7.0 nan 8.0
@@ -770,9 +770,9 @@ def merge(
     <xarray.Dataset>
     Dimensions:  (lat: 3, lon: 3, time: 2)
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0 150.0
-    * time     (time) float64 30.0 60.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0 150.0
+      * time     (time) float64 30.0 60.0
     Data variables:
         var1     (lat, lon) float64 1.0 2.0 -999.0 3.0 ... -999.0 -999.0 -999.0
         var2     (lat, lon) float64 5.0 -999.0 6.0 -999.0 ... -999.0 7.0 -999.0 8.0
@@ -782,9 +782,9 @@ def merge(
     <xarray.Dataset>
     Dimensions:  (lat: 2, lon: 2, time: 2)
     Coordinates:
-    * lat      (lat) float64 35.0 40.0
-    * lon      (lon) float64 100.0 120.0
-    * time     (time) float64 30.0 60.0
+      * lat      (lat) float64 35.0 40.0
+      * lon      (lon) float64 100.0 120.0
+      * time     (time) float64 30.0 60.0
     Data variables:
         var1     (lat, lon) float64 1.0 2.0 3.0 5.0
         var2     (lat, lon) float64 5.0 6.0 7.0 8.0
@@ -794,9 +794,9 @@ def merge(
     <xarray.Dataset>
     Dimensions:  (lat: 1, lon: 1, time: 2)
     Coordinates:
-    * lat      (lat) float64 35.0
-    * lon      (lon) float64 100.0
-    * time     (time) float64 30.0 60.0
+      * lat      (lat) float64 35.0
+      * lon      (lon) float64 100.0
+      * time     (time) float64 30.0 60.0
     Data variables:
         var1     (lat, lon) float64 1.0
         var2     (lat, lon) float64 5.0
@@ -806,9 +806,9 @@ def merge(
     <xarray.Dataset>
     Dimensions:  (lat: 1, lon: 1, time: 2)
     Coordinates:
-    * lat      (lat) float64 35.0
-    * lon      (lon) float64 100.0
-    * time     (time) float64 30.0 60.0
+      * lat      (lat) float64 35.0
+      * lon      (lon) float64 100.0
+      * time     (time) float64 30.0 60.0
     Data variables:
         var1     (lat, lon) float64 1.0
         var2     (lat, lon) float64 5.0
@@ -818,9 +818,9 @@ def merge(
     <xarray.Dataset>
     Dimensions:  (lat: 3, lon: 3, time: 2)
     Coordinates:
-    * lat      (lat) float64 35.0 40.0 42.0
-    * lon      (lon) float64 100.0 120.0 150.0
-    * time     (time) float64 30.0 60.0
+      * lat      (lat) float64 35.0 40.0 42.0
+      * lon      (lon) float64 100.0 120.0 150.0
+      * time     (time) float64 30.0 60.0
     Data variables:
         var1     (lat, lon) float64 1.0 2.0 nan 3.0 5.0 nan nan nan nan
         var2     (lat, lon) float64 5.0 nan 6.0 nan nan nan 7.0 nan 8.0
diff --git a/xarray/core/nputils.py b/xarray/core/nputils.py
index bae94121d42..c65c22f5384 100644
--- a/xarray/core/nputils.py
+++ b/xarray/core/nputils.py
@@ -173,14 +173,19 @@ def _rolling_window(a, window, axis=-1):
     Examples
     --------
     >>> x = np.arange(10).reshape((2, 5))
-    >>> np.rolling_window(x, 3, axis=-1)
-    array([[[0, 1, 2], [1, 2, 3], [2, 3, 4]],
-           [[5, 6, 7], [6, 7, 8], [7, 8, 9]]])
+    >>> _rolling_window(x, 3, axis=-1)
+    array([[[0, 1, 2],
+            [1, 2, 3],
+            [2, 3, 4]],
+    <BLANKLINE>
+           [[5, 6, 7],
+            [6, 7, 8],
+            [7, 8, 9]]])
 
     Calculate rolling mean of last dimension:
-    >>> np.mean(np.rolling_window(x, 3, axis=-1), -1)
-    array([[ 1.,  2.,  3.],
-           [ 6.,  7.,  8.]])
+    >>> np.mean(_rolling_window(x, 3, axis=-1), -1)
+    array([[1., 2., 3.],
+           [6., 7., 8.]])
 
     This function is taken from https://github.com/numpy/numpy/pull/31
     but slightly modified to accept axis option.
diff --git a/xarray/core/options.py b/xarray/core/options.py
index bb1b1c47840..5a78aa10b90 100644
--- a/xarray/core/options.py
+++ b/xarray/core/options.py
@@ -114,13 +114,14 @@ class set_options:
     <xarray.Dataset>
     Dimensions:  (x: 1000)
     Coordinates:
-      * x        (x) int64 0 1 2 3 4 5 6 ...
+      * x        (x) int64 0 1 2 ... 998 999
     Data variables:
         *empty*
 
     Or to set global options:
 
-    >>> xr.set_options(display_width=80)
+    >>> xr.set_options(display_width=80)  # doctest: +ELLIPSIS
+    <xarray.core.options.set_options object at 0x...>
     """
 
     def __init__(self, **kwargs):
diff --git a/xarray/core/rolling.py b/xarray/core/rolling.py
index 0c4614e0b57..ea99529c078 100644
--- a/xarray/core/rolling.py
+++ b/xarray/core/rolling.py
@@ -250,15 +250,29 @@ def construct(
         >>> rolling = da.rolling(b=3)
         >>> rolling.construct("window_dim")
         <xarray.DataArray (a: 2, b: 4, window_dim: 3)>
-        array([[[np.nan, np.nan, 0], [np.nan, 0, 1], [0, 1, 2], [1, 2, 3]],
-               [[np.nan, np.nan, 4], [np.nan, 4, 5], [4, 5, 6], [5, 6, 7]]])
+        array([[[nan, nan,  0.],
+                [nan,  0.,  1.],
+                [ 0.,  1.,  2.],
+                [ 1.,  2.,  3.]],
+        <BLANKLINE>
+               [[nan, nan,  4.],
+                [nan,  4.,  5.],
+                [ 4.,  5.,  6.],
+                [ 5.,  6.,  7.]]])
         Dimensions without coordinates: a, b, window_dim
 
         >>> rolling = da.rolling(b=3, center=True)
         >>> rolling.construct("window_dim")
         <xarray.DataArray (a: 2, b: 4, window_dim: 3)>
-        array([[[np.nan, 0, 1], [0, 1, 2], [1, 2, 3], [2, 3, np.nan]],
-               [[np.nan, 4, 5], [4, 5, 6], [5, 6, 7], [6, 7, np.nan]]])
+        array([[[nan,  0.,  1.],
+                [ 0.,  1.,  2.],
+                [ 1.,  2.,  3.],
+                [ 2.,  3., nan]],
+        <BLANKLINE>
+               [[nan,  4.,  5.],
+                [ 4.,  5.,  6.],
+                [ 5.,  6.,  7.],
+                [ 6.,  7., nan]]])
         Dimensions without coordinates: a, b, window_dim
 
         """
@@ -311,8 +325,15 @@ def reduce(self, func, **kwargs):
         >>> rolling = da.rolling(b=3)
         >>> rolling.construct("window_dim")
         <xarray.DataArray (a: 2, b: 4, window_dim: 3)>
-        array([[[np.nan, np.nan, 0], [np.nan, 0, 1], [0, 1, 2], [1, 2, 3]],
-               [[np.nan, np.nan, 4], [np.nan, 4, 5], [4, 5, 6], [5, 6, 7]]])
+        array([[[nan, nan,  0.],
+                [nan,  0.,  1.],
+                [ 0.,  1.,  2.],
+                [ 1.,  2.,  3.]],
+        <BLANKLINE>
+               [[nan, nan,  4.],
+                [nan,  4.,  5.],
+                [ 4.,  5.,  6.],
+                [ 5.,  6.,  7.]]])
         Dimensions without coordinates: a, b, window_dim
 
         >>> rolling.reduce(np.sum)
@@ -326,7 +347,7 @@ def reduce(self, func, **kwargs):
         <xarray.DataArray (a: 2, b: 4)>
         array([[ 0.,  1.,  3.,  6.],
                [ 4.,  9., 15., 18.]])
-
+        Dimensions without coordinates: a, b
         """
         rolling_dim = {
             d: utils.get_temp_dimname(self.obj.dims, f"_rolling_dim_{d}")
diff --git a/xarray/core/rolling_exp.py b/xarray/core/rolling_exp.py
index 96444f0f864..31d3238e978 100644
--- a/xarray/core/rolling_exp.py
+++ b/xarray/core/rolling_exp.py
@@ -93,7 +93,7 @@ def mean(self):
         >>> da = xr.DataArray([1, 1, 2, 2, 2], dims="x")
         >>> da.rolling_exp(x=2, window_type="span").mean()
         <xarray.DataArray (x: 5)>
-        array([1.      , 1.      , 1.692308, 1.9     , 1.966942])
+        array([1.        , 1.        , 1.69230769, 1.9       , 1.96694215])
         Dimensions without coordinates: x
         """
 
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index ce6df5282c5..b3e67b2f5dc 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -920,7 +920,7 @@ def copy(self, deep=True, data=None):
 
         >>> var.copy(data=[0.1, 0.2, 0.3])
         <xarray.Variable (x: 3)>
-        array([ 0.1,  0.2,  0.3])
+        array([0.1, 0.2, 0.3])
         >>> var
         <xarray.Variable (x: 3)>
         array([7, 2, 3])
@@ -1958,15 +1958,29 @@ def rolling_window(
         Examples
         --------
         >>> v = Variable(("a", "b"), np.arange(8).reshape((2, 4)))
-        >>> v.rolling_window(x, "b", 3, "window_dim")
+        >>> v.rolling_window("b", 3, "window_dim")
         <xarray.Variable (a: 2, b: 4, window_dim: 3)>
-        array([[[nan, nan, 0], [nan, 0, 1], [0, 1, 2], [1, 2, 3]],
-               [[nan, nan, 4], [nan, 4, 5], [4, 5, 6], [5, 6, 7]]])
-
-        >>> v.rolling_window(x, "b", 3, "window_dim", center=True)
+        array([[[nan, nan,  0.],
+                [nan,  0.,  1.],
+                [ 0.,  1.,  2.],
+                [ 1.,  2.,  3.]],
+        <BLANKLINE>
+               [[nan, nan,  4.],
+                [nan,  4.,  5.],
+                [ 4.,  5.,  6.],
+                [ 5.,  6.,  7.]]])
+
+        >>> v.rolling_window("b", 3, "window_dim", center=True)
         <xarray.Variable (a: 2, b: 4, window_dim: 3)>
-        array([[[nan, 0, 1], [0, 1, 2], [1, 2, 3], [2, 3, nan]],
-               [[nan, 4, 5], [4, 5, 6], [5, 6, 7], [6, 7, nan]]])
+        array([[[nan,  0.,  1.],
+                [ 0.,  1.,  2.],
+                [ 1.,  2.,  3.],
+                [ 2.,  3., nan]],
+        <BLANKLINE>
+               [[nan,  4.,  5.],
+                [ 4.,  5.,  6.],
+                [ 5.,  6.,  7.],
+                [ 6.,  7., nan]]])
         """
         if fill_value is dtypes.NA:  # np.nan is passed
             dtype, fill_value = dtypes.maybe_promote(self.dtype)