diff --git a/docs/conf.py b/docs/conf.py
index d7c8da5898..e9ffddb05c 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -2,7 +2,6 @@
 
 from __future__ import annotations
 
-import os
 import shutil
 import sys
 from datetime import datetime
@@ -17,15 +16,11 @@
 from packaging.version import Version
 from sphinxcontrib.katex import NODEJS_BINARY
 
-# Don’t use tkinter agg when importing scanpy → … → matplotlib
-matplotlib.use("agg")
+if TYPE_CHECKING:
+    from sphinx.application import Sphinx
 
 HERE = Path(__file__).parent
 sys.path[:0] = [str(HERE.parent), str(HERE / "extensions")]
-os.environ["SPHINX_RUNNING"] = "1"  # for scanpy._singleton
-
-if TYPE_CHECKING:
-    from sphinx.application import Sphinx
 
 
 # -- General configuration ------------------------------------------------
@@ -73,9 +68,9 @@
     "sphinx.ext.coverage",
     "sphinx.ext.napoleon",
     "sphinx.ext.autosummary",
+    "sphinx_exec_jupyter",
     "sphinxcontrib.bibtex",
     "sphinxcontrib.katex",
-    "matplotlib.sphinxext.plot_directive",
     "sphinx_autodoc_typehints",  # needs to be after napoleon
     "git_ref",  # needs to be before scanpydoc.rtd_github_links
     "scanpydoc",  # needs to be before sphinx.ext.linkcode
@@ -105,6 +100,14 @@
 napoleon_custom_sections = [("Params", "Parameters")]
 todo_include_todos = False
 api_dir = HERE / "api"  # function_images
+exec_jupyter_code = """
+# setup notebook backend
+import matplotlib
+matplotlib.use("module://matplotlib_inline.backend_inline")
+# import all slow optional imports before running code
+import scanpy, umap, seaborn, sklearn.metrics, pynndescent, networkx
+del scanpy, umap, seaborn, sklearn, pynndescent, networkx, matplotlib
+"""
 myst_enable_extensions = [
     "amsmath",
     "colon_fence",
@@ -119,10 +122,12 @@
     "application/vnd.microsoft.datawrangler.viewer.v0+json",
 ]
 nb_output_stderr = "remove"
-nb_execution_mode = "off"
+nb_execution_mode = "cache"
+nb_execution_excludepatterns = [
+    f"{d}{'/*' * n}" for d in ["tutorials", "how-to"] for n in (1, 2, 3)
+]
 nb_merge_streams = True
 
-
 ogp_site_url = "https://scanpy.scverse.org/en/stable/"
 ogp_image = f"{ogp_site_url}_static/Scanpy_Logo_BrightFG.svg"
 
diff --git a/hatch.toml b/hatch.toml
index 564fe0c0be..febdf1a48c 100644
--- a/hatch.toml
+++ b/hatch.toml
@@ -3,7 +3,7 @@ installer = "uv"
 dependency-groups = [ "dev" ]
 
 [envs.docs]
-dependency-groups = [ "doc" ]
+dependency-groups = [ "docs" ]
 extra-dependencies = [ "pandas>=3" ]
 scripts.build = "sphinx-build -M html docs docs/_build -W {args}"
 scripts.open = "python3 -m webbrowser -t docs/_build/html/index.html"
diff --git a/pyproject.toml b/pyproject.toml
index 16cf136060..da05accbf4 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -115,7 +115,7 @@ test = [
     "zarr>=2.18.7",
     { include-group = "test-min" },
 ]
-doc = [
+docs = [
     "ipython>=8.27",                        # for nbsphinx code highlighting
     "myst-nb>=1.4",
     "myst-parser>=2",
@@ -131,6 +131,7 @@ doc = [
     "sphinx-book-theme>=1.1",
     "sphinx-copybutton",
     "sphinx-design",
+    "sphinx-exec-jupyter>=0.2.1",
     "sphinx-issues>=5.0.1",
     "sphinxcontrib-bibtex",
     "sphinxcontrib-katex",
diff --git a/src/scanpy/metrics/_metrics.py b/src/scanpy/metrics/_metrics.py
index c3a9f86653..29c30d36df 100644
--- a/src/scanpy/metrics/_metrics.py
+++ b/src/scanpy/metrics/_metrics.py
@@ -53,7 +53,7 @@ def confusion_matrix(
     Examples
     --------
 
-    .. plot::
+    ..  exec-jupyter::
 
         import scanpy as sc; import seaborn as sns
         pbmc = sc.datasets.pbmc68k_reduced()
diff --git a/src/scanpy/plotting/_anndata.py b/src/scanpy/plotting/_anndata.py
index cc12d51562..5ea8e9614b 100755
--- a/src/scanpy/plotting/_anndata.py
+++ b/src/scanpy/plotting/_anndata.py
@@ -161,8 +161,7 @@ def scatter(  # noqa: PLR0913
     --------
     Plot two `.obs` annotations against each other and colour by a third.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -170,15 +169,13 @@ def scatter(  # noqa: PLR0913
 
     Plot expression of two genes against each other.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.scatter(adata, x="CD79A", y="CD3D", color="bulk_labels")
 
     Use a precomputed embedding via the `basis` argument.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.scatter(adata, basis="umap", color="bulk_labels")
 
@@ -618,8 +615,7 @@ def ranking(  # noqa: PLR0912, PLR0913
     PCA in :func:`~scanpy.datasets.pbmc68k_reduced` was computed on highly-variable
     genes only, so we subset to those genes before ranking.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -628,8 +624,7 @@ def ranking(  # noqa: PLR0912, PLR0913
 
     Include the lowest-loading genes alongside the highest.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.ranking(adata_hv, attr="varm", keys="PCs", indices=[0, 1, 2], include_lowest=True)
 
@@ -801,8 +796,7 @@ def violin(  # noqa: PLR0912, PLR0913, PLR0915
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -810,15 +804,13 @@ def violin(  # noqa: PLR0912, PLR0913, PLR0915
 
     Plot by category. Rotate x-axis labels so that they do not overlap.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.violin(adata, keys='S_score', groupby='bulk_labels', rotation=90)
 
     Set order of categories to be plotted or select specific categories to be plotted.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         groupby_order = ['CD34+', 'CD19+ B']
         sc.pl.violin(adata, keys='S_score', groupby='bulk_labels', rotation=90,
@@ -826,16 +818,14 @@ def violin(  # noqa: PLR0912, PLR0913, PLR0915
 
     Plot multiple keys.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.violin(adata, keys=['S_score', 'G2M_score'], groupby='bulk_labels',
             rotation=90)
 
     For large datasets consider omitting the overlaid scatter plot.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.violin(adata, keys='S_score', stripplot=False)
 
@@ -1022,15 +1012,13 @@ def clustermap(
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.krumsiek11()
         sc.pl.clustermap(adata)
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.clustermap(adata, obs_keys='cell_type')
 
@@ -1126,8 +1114,7 @@ def heatmap(  # noqa: PLR0912, PLR0913, PLR0915
 
     Examples
     --------
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -1495,8 +1482,7 @@ def tracksplot(  # noqa: PLR0912, PLR0913, PLR0915
     --------
     Using var_names as list:
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -1505,8 +1491,7 @@ def tracksplot(  # noqa: PLR0912, PLR0913, PLR0915
 
     Using var_names as dict:
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         markers = {{'T-cell': 'CD3D', 'B-cell': 'CD79A', 'myeloid': 'CST3'}}
         sc.pl.tracksplot(adata, markers, groupby='bulk_labels', dendrogram=True)
@@ -1754,8 +1739,7 @@ def dendrogram(
 
     Examples
     --------
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -1829,8 +1813,7 @@ def correlation_matrix(  # noqa: PLR0912, PLR0913, PLR0915
     --------
     Plot correlation matrix between cell type groups.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
diff --git a/src/scanpy/plotting/_dotplot.py b/src/scanpy/plotting/_dotplot.py
index 1ba66aea4e..9135b2c5e6 100644
--- a/src/scanpy/plotting/_dotplot.py
+++ b/src/scanpy/plotting/_dotplot.py
@@ -1063,8 +1063,7 @@ def dotplot(  # noqa: PLR0913
     Create a dot plot using the given markers and the PBMC example dataset grouped by
     the category `'bulk_labels'`.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -1073,8 +1072,7 @@ def dotplot(  # noqa: PLR0913
 
     Grouping `var_names` as well and specifying group colors for `groupby`:
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         from matplotlib import cm
         markers = {{'T-cell': 'CD3D', 'B-cell': 'CD79A', 'myeloid': 'CST3'}}
@@ -1083,8 +1081,7 @@ def dotplot(  # noqa: PLR0913
 
     Get `DotPlot` object for fine tuning
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         dp = sc.pl.dotplot(adata, markers, 'bulk_labels', return_fig=True)
         dp.add_totals().style(dot_edge_color='black', dot_edge_lw=0.5).show()
diff --git a/src/scanpy/plotting/_matrixplot.py b/src/scanpy/plotting/_matrixplot.py
index 70dc9fd831..00b733b6d2 100644
--- a/src/scanpy/plotting/_matrixplot.py
+++ b/src/scanpy/plotting/_matrixplot.py
@@ -58,8 +58,10 @@ class MatrixPlot(BasePlot):
 
     See Also
     --------
-    :func:`~scanpy.pl.matrixplot`: Simpler way to call MatrixPlot but with less options.
-    :func:`~scanpy.pl.rank_genes_groups_matrixplot`: to plot marker genes identified
+    :func:`~scanpy.pl.matrixplot`
+        Simpler way to call MatrixPlot but with less options.
+    :func:`~scanpy.pl.rank_genes_groups_matrixplot`
+        to plot marker genes identified
         using the :func:`~scanpy.tl.rank_genes_groups` function.
 
     Examples
@@ -67,8 +69,7 @@ class MatrixPlot(BasePlot):
     Simple visualization of the average expression of a few genes grouped by
     the category 'bulk_labels'.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -78,8 +79,7 @@ class MatrixPlot(BasePlot):
     Same visualization but passing var_names as dict, which adds a grouping of
     the genes on top of the image:
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         markers = {{'T-cell': 'CD3D', 'B-cell': 'CD79A', 'myeloid': 'CST3'}}
         sc.pl.MatrixPlot(adata, markers, groupby='bulk_labels').show()
@@ -201,8 +201,7 @@ def style(
         Examples
         --------
 
-        .. plot::
-            :context: close-figs
+        ..  exec-jupyter::
 
             import scanpy as sc
 
@@ -212,8 +211,7 @@ def style(
         Change color map and turn off edges:
 
 
-        .. plot::
-            :context: close-figs
+        ..  exec-jupyter::
 
             (
                 sc.pl.MatrixPlot(adata, markers, groupby='bulk_labels')
@@ -358,8 +356,7 @@ def matrixplot(  # noqa: PLR0913
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -368,24 +365,21 @@ def matrixplot(  # noqa: PLR0913
 
     Using var_names as dict:
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         markers = {{'T-cell': 'CD3D', 'B-cell': 'CD79A', 'myeloid': 'CST3'}}
         sc.pl.matrixplot(adata, markers, groupby='bulk_labels', dendrogram=True)
 
     Get Matrix object for fine tuning:
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         mp = sc.pl.matrixplot(adata, markers, 'bulk_labels', return_fig=True)
         mp.add_totals().style(edge_color='black').show()
 
     The axes used can be obtained using the get_axes() method
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         axes_dict = mp.get_axes()
 
diff --git a/src/scanpy/plotting/_preprocessing.py b/src/scanpy/plotting/_preprocessing.py
index db4c2afe45..db9a1eef19 100644
--- a/src/scanpy/plotting/_preprocessing.py
+++ b/src/scanpy/plotting/_preprocessing.py
@@ -45,8 +45,7 @@ def highly_variable_genes(  # noqa: PLR0912
     --------
     Compute and plot highly variable genes from raw PBMC data.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc3k()
@@ -57,8 +56,7 @@ def highly_variable_genes(  # noqa: PLR0912
 
     Plot on logarithmic axes.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.highly_variable_genes(adata, log=True)
 
diff --git a/src/scanpy/plotting/_qc.py b/src/scanpy/plotting/_qc.py
index 4c202e92a8..535f871298 100644
--- a/src/scanpy/plotting/_qc.py
+++ b/src/scanpy/plotting/_qc.py
@@ -71,8 +71,7 @@ def highest_expr_genes(
 
     Examples
     --------
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc3k()
@@ -80,8 +79,7 @@ def highest_expr_genes(
 
     Show only the top 10 genes
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.highest_expr_genes(adata, n_top=10)
 
diff --git a/src/scanpy/plotting/_stacked_violin.py b/src/scanpy/plotting/_stacked_violin.py
index a6f2ee6b59..4f764ba8ff 100644
--- a/src/scanpy/plotting/_stacked_violin.py
+++ b/src/scanpy/plotting/_stacked_violin.py
@@ -708,8 +708,7 @@ def stacked_violin(  # noqa: PLR0913
     --------
     Visualization of violin plots of a few genes grouped by the category `bulk_labels`:
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -719,16 +718,14 @@ def stacked_violin(  # noqa: PLR0913
     Same visualization but passing var_names as dict, which adds a grouping of
     the genes on top of the image:
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         markers = {{'T-cell': 'CD3D', 'B-cell': 'CD79A', 'myeloid': 'CST3'}}
         sc.pl.stacked_violin(adata, markers, groupby='bulk_labels', dendrogram=True)
 
     Get StackedViolin object for fine tuning
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         vp = sc.pl.stacked_violin(adata, markers, 'bulk_labels', return_fig=True)
         vp.add_totals().style(ylim=(0,5)).show()
diff --git a/src/scanpy/plotting/_tools/__init__.py b/src/scanpy/plotting/_tools/__init__.py
index e41f4ea908..aff7bde1c4 100644
--- a/src/scanpy/plotting/_tools/__init__.py
+++ b/src/scanpy/plotting/_tools/__init__.py
@@ -81,8 +81,7 @@ def pca_overview(adata: AnnData, **params):
 
     Examples
     --------
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc3k_processed()
@@ -132,16 +131,14 @@ def pca_loadings(
 
     Examples
     --------
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc3k_processed()
 
     Show first 3 components loadings
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.pca_loadings(adata, components = '1,2,3')
 
@@ -202,8 +199,7 @@ def pca_variance_ratio(
     --------
     Plot the variance ratio for the first 30 PCs.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc3k_processed()
@@ -211,8 +207,7 @@ def pca_variance_ratio(
 
     Plot on a logarithmic scale.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.pca_variance_ratio(adata, log=True)
 
@@ -382,8 +377,7 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -392,8 +386,7 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
 
     Plot top 10 genes (default 20 genes)
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups(adata, n_genes=10)
 
@@ -713,8 +706,7 @@ def rank_genes_groups_heatmap(
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -723,15 +715,13 @@ def rank_genes_groups_heatmap(
 
     Show gene names per group on the heatmap
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups_heatmap(adata, show_gene_labels=True)
 
     Plot top 5 genes per group (default 10 genes)
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups_heatmap(adata, n_genes=5, show_gene_labels=True)
 
@@ -786,8 +776,7 @@ def rank_genes_groups_tracksplot(
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -861,8 +850,7 @@ def rank_genes_groups_dotplot(  # noqa: PLR0913
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -870,15 +858,13 @@ def rank_genes_groups_dotplot(  # noqa: PLR0913
 
     Plot top 2 genes per group.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups_dotplot(adata,n_genes=2)
 
     Plot with scaled expressions for easier identification of differences.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups_dotplot(adata, n_genes=2, standard_scale='var')
 
@@ -887,8 +873,7 @@ def rank_genes_groups_dotplot(  # noqa: PLR0913
     and maximum values to plot are set to -4 and 4 respectively.
     Also, only genes with a log fold change of 3 or more are shown.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups_dotplot(
             adata,
@@ -903,8 +888,7 @@ def rank_genes_groups_dotplot(  # noqa: PLR0913
     Also, the last genes can be plotted. This can be useful to identify genes
     that are lowly expressed in a group. For this `n_genes=-4` is used
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups_dotplot(
             adata,
@@ -920,8 +904,7 @@ def rank_genes_groups_dotplot(  # noqa: PLR0913
     A list specific genes can be given to check their log fold change. If a
     dictionary, the dictionary keys will be added as labels in the plot.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         var_names = {{'T-cell': ['CD3D', 'CD3E', 'IL32'],
                       'B-cell': ['CD79A', 'CD79B', 'MS4A1'],
@@ -1001,8 +984,7 @@ def rank_genes_groups_stacked_violin(  # noqa: PLR0913
     --------
     Plot top marker genes per group as a stacked violin.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -1077,8 +1059,7 @@ def rank_genes_groups_matrixplot(  # noqa: PLR0913
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -1090,8 +1071,7 @@ def rank_genes_groups_matrixplot(  # noqa: PLR0913
     Also, only genes with a log fold change of 3 or more are shown.
 
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups_matrixplot(
             adata,
@@ -1107,8 +1087,7 @@ def rank_genes_groups_matrixplot(  # noqa: PLR0913
     Also, the last genes can be plotted. This can be useful to identify genes
     that are lowly expressed in a group. For this `n_genes=-4` is used
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.rank_genes_groups_matrixplot(
             adata,
@@ -1124,8 +1103,7 @@ def rank_genes_groups_matrixplot(  # noqa: PLR0913
     A list specific genes can be given to check their log fold change. If a
     dictionary, the dictionary keys will be added as labels in the plot.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         var_names = {{"T-cell": ['CD3D', 'CD3E', 'IL32'],
                       'B-cell': ['CD79A', 'CD79B', 'MS4A1'],
@@ -1216,8 +1194,7 @@ def rank_genes_groups_violin(  # noqa: PLR0913
     --------
     Plot violin distributions of top-ranked genes per group.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -1444,8 +1421,7 @@ def embedding_density(  # noqa: PLR0912, PLR0913, PLR0915
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -1454,15 +1430,13 @@ def embedding_density(  # noqa: PLR0912, PLR0913, PLR0915
 
     Plot all categories be default
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.embedding_density(adata, basis='umap', key='umap_density_phase')
 
     Plot selected categories
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.embedding_density(
             adata,
diff --git a/src/scanpy/plotting/_tools/paga.py b/src/scanpy/plotting/_tools/paga.py
index d5100f1338..2729fa1992 100644
--- a/src/scanpy/plotting/_tools/paga.py
+++ b/src/scanpy/plotting/_tools/paga.py
@@ -101,8 +101,7 @@ def paga_compare(  # noqa: PLR0912, PLR0913
     --------
     Compute a PAGA graph on the bundled PBMC dataset and show it next to the UMAP embedding.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -516,8 +515,7 @@ def paga(  # noqa: PLR0912, PLR0913, PLR0915
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc3k_processed()
@@ -526,8 +524,7 @@ def paga(  # noqa: PLR0912, PLR0913, PLR0915
 
     You can increase node and edge sizes by specifying additional arguments.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.paga(adata, node_size_scale=10, edge_width_scale=2)
 
diff --git a/src/scanpy/plotting/_tools/scatterplots.py b/src/scanpy/plotting/_tools/scatterplots.py
index e18305140f..1c57415fb1 100644
--- a/src/scanpy/plotting/_tools/scatterplots.py
+++ b/src/scanpy/plotting/_tools/scatterplots.py
@@ -124,8 +124,7 @@ def embedding(  # noqa: PLR0912, PLR0913, PLR0915
     --------
     Plot a precomputed UMAP embedding coloured by the `'bulk_labels'` cell-type annotation.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -134,8 +133,7 @@ def embedding(  # noqa: PLR0912, PLR0913, PLR0915
     Show several panels in a single call by passing a list of keys to `color`,
     mixing categorical annotations and gene expression.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.embedding(adata, basis="umap", color=["bulk_labels", "CD3D", "LYZ"])
 
@@ -689,8 +687,7 @@ def umap(adata: AnnData, **kwargs) -> Figure | Axes | list[Axes] | None:
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -698,22 +695,19 @@ def umap(adata: AnnData, **kwargs) -> Figure | Axes | list[Axes] | None:
 
     Colour points by discrete variable (Louvain clusters).
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.umap(adata, color="louvain")
 
     Colour points by gene expression.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.umap(adata, color="HES4")
 
     Plot muliple umaps for different gene expressions.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.umap(adata, color=["HES4", "TNFRSF4"])
 
@@ -750,8 +744,7 @@ def tsne(adata: AnnData, **kwargs) -> Figure | Axes | list[Axes] | None:
 
     Examples
     --------
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -789,8 +782,7 @@ def diffmap(adata: AnnData, **kwargs) -> Figure | Axes | list[Axes] | None:
 
     Examples
     --------
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -835,8 +827,7 @@ def draw_graph(
 
     Examples
     --------
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -893,8 +884,7 @@ def pca(
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc3k_processed()
@@ -902,15 +892,13 @@ def pca(
 
     Colour points by discrete variable (Louvain clusters).
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.pca(adata, color="louvain")
 
     Colour points by gene expression.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.pca(adata, color="CST3")
 
diff --git a/src/scanpy/preprocessing/_qc.py b/src/scanpy/preprocessing/_qc.py
index 04340ec7e9..26e6030af6 100644
--- a/src/scanpy/preprocessing/_qc.py
+++ b/src/scanpy/preprocessing/_qc.py
@@ -252,8 +252,7 @@ def calculate_qc_metrics(
     -------
     Calculate qc metrics for visualization.
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         import seaborn as sns
@@ -268,8 +267,7 @@ def calculate_qc_metrics(
             kind="hex",
         )
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sns.histplot(pbmc.obs["pct_counts_mito"])
 
diff --git a/src/scanpy/tools/_embedding_density.py b/src/scanpy/tools/_embedding_density.py
index 37e2bcf606..1c58102b69 100644
--- a/src/scanpy/tools/_embedding_density.py
+++ b/src/scanpy/tools/_embedding_density.py
@@ -88,8 +88,7 @@ def embedding_density(  # noqa: PLR0912
     Examples
     --------
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         import scanpy as sc
         adata = sc.datasets.pbmc68k_reduced()
@@ -99,8 +98,7 @@ def embedding_density(  # noqa: PLR0912
             adata, basis='umap', key='umap_density_phase', group='G1'
         )
 
-    .. plot::
-        :context: close-figs
+    ..  exec-jupyter::
 
         sc.pl.embedding_density(
             adata, basis='umap', key='umap_density_phase', group='S'