diff --git a/.gitignore b/.gitignore
index 870125286..cd7ca61c4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -24,6 +24,8 @@ pip-log.txt
 doc/_build/
 doc/build/
 .jupyterlite.doit.db
+/doc/jupyter_execute/
+doc/source/regression/*.ipynb
 
 # Editors
 .idea
diff --git a/doc/source/_static/myst-nb.css b/doc/source/_static/myst-nb.css
new file mode 100644
index 000000000..bd30038a1
--- /dev/null
+++ b/doc/source/_static/myst-nb.css
@@ -0,0 +1,82 @@
+/* MyST-NB
+
+This stylesheet targets elements output by MyST-NB that represent notebook
+cells.
+
+In some cases these rules override MyST-NB. In some cases they override PyData
+Sphinx Theme or Sphinx. And in some cases they do not override existing styling
+but add new styling. */
+
+/* Set up a few variables for this stylesheet */
+.cell,
+.pywt-handcoded-cell-output {
+  --pywt-cell-input-border-left-width: .2rem;
+
+  /* This matches the padding applied to <pre> elements by PyData Sphinx Theme */
+  --pywt-code-block-padding: 1rem;
+
+  /* override mystnb */
+  --mystnb-source-border-radius: .25rem; /* match PyData Sphinx Theme */
+}
+
+.cell .cell_input::before {
+  content: "In";
+  border-bottom: var(--mystnb-source-border-width) solid var(--pst-color-border);
+  font-weight: var(--pst-font-weight-caption);
+
+  /* Left-aligns the text in this box and the one that follows it */
+  padding-left: var(--pywt-code-block-padding);
+}
+
+/* Cannot use `.cell .cell_input` selector because the stylesheet from MyST-NB
+   uses `div.cell div.cell_input` and we want to override those rules */
+div.cell div.cell_input {
+  background-color: inherit;
+  border-color: var(--pst-color-border);
+  border-left-width: var(--pywt-cell-input-border-left-width);
+  background-clip: padding-box;
+  overflow: hidden;
+}
+
+.cell .cell_output,
+.pywt-handcoded-cell-output {
+  border: var(--mystnb-source-border-width) solid var(--pst-color-surface);
+  border-radius: var(--mystnb-source-border-radius);
+  background-clip: padding-box;
+  overflow: hidden;
+}
+
+.cell .cell_output::before,
+.pywt-handcoded-cell-output::before {
+  content: "Out";
+  display: block;
+  font-weight: var(--pst-font-weight-caption);
+
+  /* Left-aligns the text in this box and the one that follows it */
+  padding-left: var(--pywt-code-block-padding);
+}
+
+.cell .cell_output .output {
+  background-color: inherit;
+  border: none;
+  margin-top: 0;
+}
+
+.cell .cell_output,
+/* must prefix the following selector with `div.` to override Sphinx margin rule on div[class*=highlight-] */
+div.pywt-handcoded-cell-output {
+  /* Left-align the text in the output with the text in the input */
+  margin-left: calc(var(--pywt-cell-input-border-left-width) - var(--mystnb-source-border-width));
+}
+
+.cell .cell_output .output,
+.cell .cell_input pre,
+.cell .cell_output pre,
+.pywt-handcoded-cell-output .highlight,
+.pywt-handcoded-cell-output pre {
+  border-radius: 0;
+}
+
+.pywt-handcoded-cell-output pre {
+  border: none; /* MyST-NB sets border to none for <pre> tags inside div.cell */
+}
diff --git a/doc/source/_static/pywavelets.css b/doc/source/_static/pywavelets.css
index 56f12e4df..7aa4e8326 100644
--- a/doc/source/_static/pywavelets.css
+++ b/doc/source/_static/pywavelets.css
@@ -1,27 +1,37 @@
 /* Custom CSS rules for the interactive documentation button */
 
 .try_examples_button {
-    color: white;
-    background-color: #0054a6;
-    border: none;
-    padding: 5px 10px;
-    border-radius: 10px;
-    margin-bottom: 5px;
-    box-shadow: 0 2px 5px rgba(108,108,108,0.2);
-    font-weight: bold;
-    font-size: small;
+  color: white;
+  background-color: #0054a6;
+  border: none;
+  padding: 5px 10px;
+  border-radius: 10px;
+  margin-bottom: 5px;
+  box-shadow: 0 2px 5px rgba(108,108,108,0.2);
+  font-weight: bold;
+  font-size: small;
 }
 
 .try_examples_button:hover {
-background-color: #0066cc;
-transform: scale(1.02);
-box-shadow: 0 2px 5px rgba(0, 0, 0, 0.2);
-cursor: pointer;
+  background-color: #0066cc;
+  transform: scale(1.02);
+  box-shadow: 0 2px 5px rgba(0, 0, 0, 0.2);
+  cursor: pointer;
 }
 
 .try_examples_button_container {
-    display: flex;
-    justify-content: flex-start;
-    gap: 10px;
-    margin-bottom: 20px;
+  display: flex;
+  justify-content: flex-start;
+  gap: 10px;
+  margin-bottom: 20px;
+}
+
+/*
+Admonitions on this site are styled with some top margin. This makes sense when
+the admonition appears within the flow of the article. But when it is the very
+first child of an article, its top margin gets added to the article's top
+padding, resulting in too much whitespace.
+*/
+.admonition.pywt-margin-top-0 {
+  margin-top: 0;
 }
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 2820d8a08..1472a694e 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -11,9 +11,10 @@
 # serve to show the default.
 
 import datetime
-import importlib.metadata
+import os
+import re
+from pathlib import Path
 
-import jinja2.filters
 import numpy as np
 
 import pywt
@@ -25,6 +26,76 @@
 except TypeError:
     pass
 
+from sphinx.application import Sphinx
+
+HERE = Path(__file__).parent
+
+
+def preprocess_notebooks(app: Sphinx, *args, **kwargs):
+    """Preprocess Markdown notebooks to convert them to IPyNB format."""
+
+    import jupytext
+    import nbformat
+
+    print("Converting Markdown files to IPyNB...")
+    for path in (HERE / "regression").glob("*.md"):
+        if any(path.match(pattern) for pattern in exclude_patterns):
+            continue
+        nb = jupytext.read(str(path))
+
+        # In .md to .ipynb conversion, do not include any cells that have the
+        # jupyterlite_sphinx_strip tag
+        nb.cells = [
+            cell for cell in nb.cells if "jupyterlite_sphinx_strip" not in cell.metadata.get("tags", [])
+        ]
+
+        ipynb_path = path.with_suffix(".ipynb")
+        with open(ipynb_path, "w") as f:
+            nbformat.write(nb, f)
+            print(f"Converted {path} to {ipynb_path}")
+
+
+# Should match {{ parent_docname }} or {{parent_docname}}
+parent_docname_substitution_re = re.compile(r"{{\s*parent_docname\s*}}")
+
+def sub_parent_docname_in_header(
+    app: Sphinx, relative_path: Path, parent_docname: str, content: list[str]
+):
+    """Fill in the name of the document in the header.
+
+    When regression/header.md is read via the include directive, replace
+    {{ parent_docname }} with the name of the parent document that included
+    header.md.
+
+    Note: parent_docname does not include the file extension.
+
+    Here is a simplified example of how this works.
+
+    Contents of header.md:
+
+        {download}`Download {{ parent_docname }}.md <{{ parent_docname }}.md>`
+
+    Contents of foobar.md:
+
+        ```{include} header.md
+        ```
+
+    After this function and others are run...
+
+    Contents of foobar.md:
+
+        {download}`Download foobar.md <foobar.md>`
+    """
+    if not relative_path.match("regression/header.md"):
+        return
+
+    for i, value in enumerate(content):
+        content[i] = re.sub(parent_docname_substitution_re, parent_docname, value)
+
+
+def setup(app):
+    app.connect("config-inited", preprocess_notebooks)
+    app.connect("include-read", sub_parent_docname_in_header)
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
@@ -38,6 +109,7 @@
 extensions = [
     'jupyterlite_sphinx',
     'matplotlib.sphinxext.plot_directive',
+    'myst_nb',
     'numpydoc',
     'sphinx.ext.autodoc',
     'sphinx.ext.autosummary',
@@ -47,15 +119,19 @@
     'sphinx.ext.mathjax',
     'sphinx.ext.todo',
     'sphinx_copybutton',
+    'sphinx_design',
     'sphinx_togglebutton',
-
 ]
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
 
 # The suffix of source filenames.
-source_suffix = '.rst'
+source_suffix = {
+    '.rst': 'restructuredtext',
+    '.md': 'myst-nb',
+    'ipynb': None,  # do not parse IPyNB files
+}
 
 # The encoding of source files.
 source_encoding = 'utf-8'
@@ -141,8 +217,8 @@
 "show_prev_next": True,
 "footer_start": ["copyright", "sphinx-version"],
 "footer_end": ["theme-version"],
-"pygment_light_style": "a11y-high-contrast-light",
-"pygment_dark_style": "a11y-high-contrast-dark",
+"pygments_light_style": "a11y-high-contrast-light",
+"pygments_dark_style": "a11y-high-contrast-dark",
 "icon_links": [
         {
             "name": "Discussion group on Google Groups",
@@ -200,6 +276,7 @@
 # _static directory.
 html_css_files = [
     "pywavelets.css",
+    "myst-nb.css"
 ]
 
 # If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
@@ -278,7 +355,12 @@
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
-exclude_patterns = ['substitutions.rst', ]
+exclude_patterns = [
+    'substitutions.rst',
+    'regression/header.md',
+    'regression/README.md',
+    'regression/*.ipynb'  # exclude IPyNB files from the build
+]
 
 # numpydoc_show_class_members = False
 numpydoc_class_members_toctree = False
@@ -310,3 +392,25 @@
 Shall you encounter any issues, please feel free to report them on the
 [PyWavelets issue tracker](https://github.com/PyWavelets/pywt/issues)."""
 )
+
+jupyterlite_silence = True
+strip_tagged_cells = True
+
+# -- Options for MyST-NB and Markdown-based content --------------------------
+
+os.environ["PYDEVD_DISABLE_FILE_VALIDATION"] = "1"
+
+# https://myst-nb.readthedocs.io/en/latest/configuration.html
+nb_execution_mode = 'auto'
+nb_execution_timeout = 60
+nb_execution_allow_errors = False
+nb_execution_raise_on_error = True
+nb_render_markdown_format = "myst"
+nb_remove_code_source = False
+nb_remove_code_outputs = False
+
+myst_enable_extensions = [
+    'amsmath',
+    'colon_fence',
+    'dollarmath',
+]
diff --git a/doc/source/regression/README.md b/doc/source/regression/README.md
new file mode 100644
index 000000000..7b1b58c96
--- /dev/null
+++ b/doc/source/regression/README.md
@@ -0,0 +1,54 @@
+# Regression folder
+
+This folder contains various useful examples illustrating how to use and how not
+to use PyWavelets.
+
+The examples are written in the [MyST markdown notebook
+format](https://myst-nb.readthedocs.io/en/v0.13.2/use/markdown.html). This
+allows each .md file to function simultaneously as documentation that can be fed
+into Sphinx and as a source file that can be converted to the Jupyter notebook
+format (.ipynb), which can then be opened in notebook applications such as
+JupyterLab. For this reason, each example page in this folder includes a header template
+that adds a blurb to the top of each page about how the page can be
+run or downloaded as a Jupyter notebook.
+
+There are a few shortcomings to this approach of generating the code cell outputs in
+the documentation pages at build time rather than hand editing them into the
+document source file. One is that we can no longer compare the generated outputs
+with the expected outputs as we used to do with doctest. Another is that we
+lose some control over how we want the outputs to appear, unless we use a workaround.
+
+Here is the workaround we created. First we tell MyST-NB to remove the generated
+cell output from the documentation page by adding the `remove-output` tag to the
+`code-cell` directive in the markdown file. Then we hand code the output in a
+`code-block` directive, not to be confused with `code-cell`! The `code-cell`
+directive says "I am notebook code cell input, run me!" The `code-block`
+directive says, "I am just a block of code for documentation purposes, don't run
+me!" To the code block, we add the `.pywt-handcoded-cell-output` class so that
+we can style it to look the same as other cell outputs on the same HTML page.
+Finally, we tag the handcoded output with `jupyterlite_sphinx_strip` so that we
+can exclude it when converting from .md to .ipynb. That way only generated
+output appears in the .ipynb notebook.
+
+To recap:
+
+- We use the `remove-output` tag to remove the **generated** code cell output
+  during .md to .html conversion (this conversion is done by MyST-NB).
+- We use the `jupyterlite_sphinx_strip` tag to remove the **handcoded** output
+  during .md to .ipynb conversion (this conversion is done by Jupytext).
+
+Example markdown:
+
+    ```{code-cell}
+    :tags: [raises-exception, remove-output]
+    1 / 0
+    ```
+
+    +++ {"tags" ["jupyterlite_sphinx_strip"]}
+
+    ```{code-block} python
+    :class: pywt-handcoded-cell-output
+    Traceback (most recent call last):
+    ...
+    ZeroDivisionError: division by zero
+    ```
diff --git a/doc/source/regression/dwt-idwt.md b/doc/source/regression/dwt-idwt.md
new file mode 100644
index 000000000..101295aca
--- /dev/null
+++ b/doc/source/regression/dwt-idwt.md
@@ -0,0 +1,212 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 1.0.0
+    jupytext_version: 1.16.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{include} header.md
+
+```
+
+# DWT and IDWT
+
+## Discrete Wavelet Transform
+
+Let's do a [Discrete Wavelet Transform](ref-dwt) of some sample data `x`
+using the `db2` wavelet. It's simple:
+
+```{code-cell}
+import pywt
+x = [3, 7, 1, 1, -2, 5, 4, 6]
+cA, cD = pywt.dwt(x, 'db2')
+```
+
+And the approximation and details coefficients are in `cA` and `cD`
+respectively:
+
+```{code-cell}
+cA
+```
+
+```{code-cell}
+cD
+```
+
+## Inverse Discrete Wavelet Transform
+
+Now let's do the opposite operation, an [Inverse Discrete Wavelet Transform](ref-idwt):
+
+```{code-cell}
+pywt.idwt(cA, cD, 'db2')
+```
+
+Voilà! That's it!
+
+## More examples
+
+Now let's experiment with `dwt()` some more. For example, let's pass a
+`Wavelet` object instead of the wavelet name and specify the signal
+extension mode (the default is `Modes.symmetric`) for the border effect
+handling:
+
+```{code-cell}
+w = pywt.Wavelet('sym3')
+cA, cD = pywt.dwt(x, wavelet=w, mode='constant')
+```
+
+```{code-cell}
+print(cA)
+```
+
+```{code-cell}
+print(cD)
+```
+
+Note that the output coefficients arrays' length depends not only on the
+input data length but also on the `Wavelet` type (particularly on its
+filters length `Wavelet.dec_len` that are used in the transformation).
+
+To find out what the size of the output data will be, use the `dwt_coeff_len()`
+function:
+
+```{code-cell}
+pywt.dwt_coeff_len(data_len=len(x), filter_len=w.dec_len, mode='symmetric')
+```
+
+```{code-cell}
+pywt.dwt_coeff_len(len(x), w, 'symmetric')
+```
+
+and the length of the output is:
+
+```{code-cell}
+len(cA)
+```
+
+Looks fine. (And if you expected that the output length would be a half of the
+input data length, well, that's the trade-off that allows for the perfect
+reconstruction...).
+
+The third argument of the `dwt_coeff_len()` function is the already mentioned signal
+extension mode (please refer to the PyWavelets' documentation for the `modes`
+description). Currently, there are six extension modes available under `Modes`:
+
+```{code-cell}
+pywt.Modes.modes
+```
+
+As you see in the above example, the periodization (`Modes.periodization`) mode
+is slightly different from the others. Its aim when doing the `pywt.dwt` transform
+is to output coefficients arrays that are half of the length of the input data.
+
+Knowing that, you should never mix the periodization mode with other modes when
+doing `dwt` and `idwt`. Otherwise, it will produce **invalid results**:
+
+```{code-cell}
+x = [3, 7, 1, 1, -2, 5, 4, 6]
+
+cA, cD = pywt.dwt(x, wavelet=w, mode='periodization')
+print(pywt.idwt(cA, cD, 'sym3', 'symmetric'))  # invalid mode
+```
+
+```{code-cell}
+print(pywt.idwt(cA, cD, 'sym3', 'periodization'))
+```
+
+## Tips & tricks
+
+### Passing `None` instead of coefficients data to `pywt.idwt()`
+
+Now, we showcase some tips and tricks. Passing `None` as one of the coefficient
+arrays parameters is similar to passing a _zero-filled_ array. The results are
+simply the same:
+
+```{code-cell}
+print(pywt.idwt([1,2,0,1], None, 'db2', 'symmetric'))
+```
+
+```{code-cell}
+print(pywt.idwt([1, 2, 0, 1], [0, 0, 0, 0], 'db2', 'symmetric'))
+```
+
+```{code-cell}
+print(pywt.idwt(None, [1, 2, 0, 1], 'db2', 'symmetric'))
+```
+
+```{code-cell}
+print(pywt.idwt([0, 0, 0, 0], [1, 2, 0, 1], 'db2', 'symmetric'))
+```
+
+Remember that only one argument at a time can be `None`:
+
+```{code-cell}
+---
+tags: [raises-exception, remove-output]
+---
+print(pywt.idwt(None, None, 'db2', 'symmetric'))
+```
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{code-block} python
+:class: pywt-handcoded-cell-output
+Traceback (most recent call last):
+...
+ValueError: At least one coefficient parameter must be specified.
+```
+
+### Coefficients data size in `pywt.idwt`
+
+When doing the `idwt` transform, usually the coefficient arrays
+must have the same size.
+
+```{code-cell}
+---
+tags: [raises-exception, remove-output]
+---
+print(pywt.idwt([1, 2, 3, 4, 5], [1, 2, 3, 4], 'db2', 'symmetric'))
+```
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{code-block} python
+:class: pywt-handcoded-cell-output
+Traceback (most recent call last):
+...
+ValueError: Coefficients arrays must have the same size.
+```
+
+Not every coefficient array can be used in `idwt`. In the
+following example the `idwt` will fail because the input arrays are
+invalid - they couldn't be created as a result of `dwt`, because
+the minimal output length for dwt using `db4` wavelet and the `Modes.symmetric`
+mode is `4`, not `3`:
+
+```{code-cell}
+---
+tags: [raises-exception, remove-output]
+---
+pywt.idwt([1,2,4], [4,1,3], 'db4', 'symmetric')
+```
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{code-block} python
+:class: pywt-handcoded-cell-output
+Traceback (most recent call last):
+...
+ValueError: Invalid coefficient arrays length for specified wavelet. Wavelet and mode must be the same as used for decomposition.
+```
+
+```{code-cell}
+int(pywt.dwt_coeff_len(1, pywt.Wavelet('db4').dec_len, 'symmetric'))
+```
diff --git a/doc/source/regression/dwt-idwt.rst b/doc/source/regression/dwt-idwt.rst
deleted file mode 100644
index a0552759b..000000000
--- a/doc/source/regression/dwt-idwt.rst
+++ /dev/null
@@ -1,151 +0,0 @@
-.. _reg-dwt-idwt:
-
-.. currentmodule:: pywt
-
-DWT and IDWT
-============
-
-Discrete Wavelet Transform
---------------------------
-
-Let's do a :func:`Discrete Wavelet Transform <dwt>` of a sample data ``x``
-using the ``db2`` wavelet. It's simple..
-
-    >>> import pywt
-    >>> x = [3, 7, 1, 1, -2, 5, 4, 6]
-    >>> cA, cD = pywt.dwt(x, 'db2')
-
-And the approximation and details coefficients are in ``cA`` and ``cD``
-respectively:
-
-    >>> print(cA)
-    [ 5.65685425  7.39923721  0.22414387  3.33677403  7.77817459]
-    >>> print(cD)
-    [-2.44948974 -1.60368225 -4.44140056 -0.41361256  1.22474487]
-
-Inverse Discrete Wavelet Transform
-----------------------------------
-
-Now let's do an opposite operation
-- :func:`Inverse Discrete Wavelet Transform <idwt>`:
-
-    >>> print(pywt.idwt(cA, cD, 'db2'))
-    [ 3.  7.  1.  1. -2.  5.  4.  6.]
-
-Voilà! That's it!
-
-More Examples
--------------
-
-Now let's experiment with the :func:`dwt` some more. For example let's pass a
-:class:`Wavelet` object instead of the wavelet name and specify signal
-extension mode (the default is :ref:`symmetric <Modes.symmetric>`) for the
-border effect handling:
-
-    >>> w = pywt.Wavelet('sym3')
-    >>> cA, cD = pywt.dwt(x, wavelet=w, mode='constant')
-    >>> print(cA)
-    [ 4.38354585  3.80302657  7.31813271 -0.58565539  4.09727044  7.81994027]
-    >>> print(cD)
-    [-1.33068221 -2.78795192 -3.16825651 -0.67715519 -0.09722957 -0.07045258]
-
-Note that the output coefficients arrays length depends not only on the input
-data length but also on the :class:Wavelet type (particularly on its
-:attr:`filters length <~Wavelet.dec_len>` that are used in the transformation).
-
-To find out what will be the output data size use the :func:`dwt_coeff_len`
-function:
-
-    >>> # int() is for normalizing Python integers and long integers for documentation tests
-    >>> int(pywt.dwt_coeff_len(data_len=len(x), filter_len=w.dec_len, mode='symmetric'))
-    6
-    >>> int(pywt.dwt_coeff_len(len(x), w, 'symmetric'))
-    6
-    >>> len(cA)
-    6
-
-Looks fine. (And if you expected that the output length would be a half of the
-input data length, well, that's the trade-off that allows for the perfect
-reconstruction...).
-
-The third argument of the :func:`dwt_coeff_len` is the already mentioned signal
-extension mode (please refer to the PyWavelets' documentation for the
-:ref:`modes <modes>` description). Currently there are six
-:ref:`extension modes <Modes>` available:
-
-    >>> pywt.Modes.modes
-    ['zero', 'constant', 'symmetric', 'periodic', 'smooth', 'periodization', 'reflect', 'antisymmetric', 'antireflect']
-
-As you see in the above example, the :ref:`periodization <Modes.periodization>`
-(periodization) mode is slightly different from the others. It's aim when
-doing the :func:`DWT <dwt>` transform is to output coefficients arrays that
-are half of the length of the input data.
-
-Knowing that, you should never mix the periodization mode with other modes when
-doing :func:`DWT <dwt>` and :func:`IDWT <idwt>`. Otherwise, it will produce
-**invalid results**:
-
-    >>> x
-    [3, 7, 1, 1, -2, 5, 4, 6]
-    >>> cA, cD = pywt.dwt(x, wavelet=w, mode='periodization')
-    >>> print(pywt.idwt(cA, cD, 'sym3', 'symmetric')) # invalid mode
-    [ 1.  1. -2.  5.]
-    >>> print(pywt.idwt(cA, cD, 'sym3', 'periodization'))
-    [ 3.  7.  1.  1. -2.  5.  4.  6.]
-
-
-Tips & tricks
--------------
-
-Passing ``None`` instead of coefficients data to :func:`idwt`
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Now some tips & tricks. Passing ``None`` as one of the coefficient arrays
-parameters is similar to passing a *zero-filled* array. The results are simply
-the same:
-
-    >>> print(pywt.idwt([1,2,0,1], None, 'db2', 'symmetric'))
-    [ 1.19006969  1.54362308  0.44828774 -0.25881905  0.48296291  0.8365163 ]
-
-    >>> print(pywt.idwt([1, 2, 0, 1], [0, 0, 0, 0], 'db2', 'symmetric'))
-    [ 1.19006969  1.54362308  0.44828774 -0.25881905  0.48296291  0.8365163 ]
-
-    >>> print(pywt.idwt(None, [1, 2, 0, 1], 'db2', 'symmetric'))
-    [ 0.57769726 -0.93125065  1.67303261 -0.96592583 -0.12940952 -0.22414387]
-
-    >>> print(pywt.idwt([0, 0, 0, 0], [1, 2, 0, 1], 'db2', 'symmetric'))
-    [ 0.57769726 -0.93125065  1.67303261 -0.96592583 -0.12940952 -0.22414387]
-
-Remember that only one argument at a time can be ``None``:
-
-    >>> print(pywt.idwt(None, None, 'db2', 'symmetric'))
-    Traceback (most recent call last):
-    ...
-    ValueError: At least one coefficient parameter must be specified.
-
-
-Coefficients data size in :attr:`idwt`
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-When doing the :func:`IDWT <idwt>` transform, usually the coefficient arrays
-must have the same size.
-
-    >>> print(pywt.idwt([1, 2, 3, 4, 5], [1, 2, 3, 4], 'db2', 'symmetric'))
-    Traceback (most recent call last):
-    ...
-    ValueError: Coefficients arrays must have the same size.
-
-
-Not every coefficient array can be used in :func:`IDWT <idwt>`. In the
-following example the :func:`idwt` will fail because the input arrays are
-invalid - they couldn't be created as a result of :func:`DWT <dwt>`, because
-the minimal output length for dwt using ``db4`` wavelet and the :ref:`symmetric
-<Modes.symmetric>` mode is ``4``, not ``3``:
-
-    >>> pywt.idwt([1,2,4], [4,1,3], 'db4', 'symmetric')
-    Traceback (most recent call last):
-    ...
-    ValueError: Invalid coefficient arrays length for specified wavelet. Wavelet and mode must be the same as used for decomposition.
-
-    >>> int(pywt.dwt_coeff_len(1, pywt.Wavelet('db4').dec_len, 'symmetric'))
-    4
diff --git a/doc/source/regression/gotchas.md b/doc/source/regression/gotchas.md
new file mode 100644
index 000000000..6cb955730
--- /dev/null
+++ b/doc/source/regression/gotchas.md
@@ -0,0 +1,53 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 1.0.0
+    jupytext_version: 1.16.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{include} header.md
+
+```
+
+# Gotchas
+
+PyWavelets utilizes `NumPy` under the hood. That's why handling the data
+that contains `None` values can be surprising. `None` values are converted to
+'not a number' (`numpy.nan`) values:
+
+```{code-cell}
+import numpy, pywt
+x = [None, None]
+mode = 'symmetric'
+wavelet = 'db1'
+cA, cD = pywt.dwt(x, wavelet, mode)
+```
+
+The results are:
+
+```{code-cell}
+numpy.all(numpy.isnan(cA))
+```
+
+<!-- True -->
+
+```{code-cell}
+numpy.all(numpy.isnan(cD))
+```
+
+<!-- True -->
+
+```{code-cell}
+rec = pywt.idwt(cA, cD, wavelet, mode)
+numpy.all(numpy.isnan(rec))
+```
+
+<!-- True -->
diff --git a/doc/source/regression/gotchas.rst b/doc/source/regression/gotchas.rst
deleted file mode 100644
index 04e267ff0..000000000
--- a/doc/source/regression/gotchas.rst
+++ /dev/null
@@ -1,25 +0,0 @@
-.. _reg-gotchas:
-
-.. currentmodule:: pywt
-
-
-=======
-Gotchas
-=======
-
-PyWavelets utilizes ``NumPy`` under the hood. That's why handling the data
-containing ``None`` values can be surprising. ``None`` values are converted to
-'not a number' (``numpy.NaN``) values:
-
-    >>> import numpy, pywt
-    >>> x = [None, None]
-    >>> mode = 'symmetric'
-    >>> wavelet = 'db1'
-    >>> cA, cD = pywt.dwt(x, wavelet, mode)
-    >>> numpy.all(numpy.isnan(cA))
-    True
-    >>> numpy.all(numpy.isnan(cD))
-    True
-    >>> rec = pywt.idwt(cA, cD, wavelet, mode)
-    >>> numpy.all(numpy.isnan(rec))
-    True
diff --git a/doc/source/regression/header.md b/doc/source/regression/header.md
new file mode 100644
index 000000000..05b713223
--- /dev/null
+++ b/doc/source/regression/header.md
@@ -0,0 +1,15 @@
+````{tip}
+:class: pywt-margin-top-0
+
+This page can also be run or downloaded as a Jupyter notebook.
+
+- Run in a new tab using JupyterLite[^what-is-lite]:
+  ```{notebooklite} {{ parent_docname }}.ipynb
+  :new_tab: True
+  ```
+- Download:
+  - Jupyter notebook: {download}`{{ parent_docname }}.ipynb <{{ parent_docname }}.ipynb>`
+  - MyST markdown: {download}`{{ parent_docname }}.md <{{ parent_docname }}.md>`
+
+[^what-is-lite]: [JupyterLite](https://jupyterlite.readthedocs.io) is a distribution of JupyterLab that runs entirely within a web browser. All computations run in the browser; none are sent to a server.
+````
diff --git a/doc/source/regression/index.rst b/doc/source/regression/index.rst
index 21065886a..1e5ab54dc 100644
--- a/doc/source/regression/index.rst
+++ b/doc/source/regression/index.rst
@@ -5,11 +5,8 @@
 Usage examples
 ==============
 
-The following examples are used as doctest regression tests written using reST
-markup. They are included in the documentation since they contain various useful
-examples illustrating how to use and how not to use PyWavelets.
-
-For more usage examples see the `demo`_ directory in the source package.
+The following pages contain various useful examples illustrating how to use and
+how not to use PyWavelets.
 
 .. toctree::
    :maxdepth: 1
@@ -22,4 +19,19 @@ For more usage examples see the `demo`_ directory in the source package.
    wp2d
    gotchas
 
+The examples in this section are written in the `MyST markdown notebook format`_. This allows
+each .md file to function simultaneously as documentation that can be fed into
+Sphinx and as a source file that can be converted to the Jupyter notebook format
+(.ipynb), which can then be opened in notebook applications such as JupyterLab.
+For this reason, each example page in this folder includes a header template
+that adds a blurb to the top of each page about how the page can be run or
+downloaded as a Jupyter notebook. This is explained in greater detail in the
+`regression folder README`_.
+
+For more usage examples see the `demo`_ directory in the source package.
+
+
 .. include:: ../common_refs.rst
+
+.. _MyST markdown notebook format: https://myst-nb.readthedocs.io/en/v0.13.2/use/markdown.html
+.. _regression folder readme: https://github.com/PyWavelets/pywt/blob/main/doc/source/regression/README.md
diff --git a/doc/source/regression/modes.md b/doc/source/regression/modes.md
new file mode 100644
index 000000000..4905593a0
--- /dev/null
+++ b/doc/source/regression/modes.md
@@ -0,0 +1,101 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 1.0.0
+    jupytext_version: 1.16.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{include} header.md
+
+```
+
+# Signal Extension Modes
+
+Import `pywt` first:
+
+```{code-cell}
+import pywt
+```
+
+and construct a helper function that can format arrays in a consistent manner
+across different systems. Please note that this function is just for the purpose of
+this example and is not part of the PyWavelets library, and it is not necessary or
+required to use it in your own code:
+
+```{code-cell}
+def format_array(a):
+    """Consistent array representation across different systems"""
+    import numpy
+    a = numpy.where(numpy.abs(a) < 1e-5, 0, a)
+    return numpy.array2string(a, precision=5, separator=' ', suppress_small=True)
+```
+
+A list of available signal extension [modes](Modes):
+
+```{code-cell}
+pywt.Modes.modes
+```
+
+An invalid mode name should raise a `ValueError`:
+
+```{code-cell}
+---
+tags: [raises-exception, remove-output]
+---
+pywt.dwt([1,2,3,4], 'db2', 'invalid')
+```
+
+```{code-block} python
+:class: pywt-handcoded-cell-output
+Traceback (most recent call last):
+...
+ValueError: Unknown mode name 'invalid'.
+```
+
+You can also refer to modes via the attributes of the `Modes` class:
+
+```{code-cell}
+x = [1, 2, 1, 5, -1, 8, 4, 6]
+for mode_name in ['zero', 'constant', 'symmetric', 'reflect', 'periodic', 'smooth', 'periodization']:
+    mode = getattr(pywt.Modes, mode_name)
+    cA, cD = pywt.dwt(x, 'db2', mode)
+    print("Mode: %d (%s)" % (mode, mode_name))
+```
+
+The default mode is [symmetric](Modes.symmetric):
+
+```{code-cell}
+cA, cD = pywt.dwt(x, 'db2')
+cA
+```
+
+```{code-cell}
+cD
+```
+
+```{code-cell}
+pywt.idwt(cA, cD, 'db2')
+```
+
+Specify the mode using a keyword argument:
+
+```{code-cell}
+cA, cD = pywt.dwt(x, 'db2', mode='symmetric')
+cA
+```
+
+```{code-cell}
+cD
+```
+
+```{code-cell}
+pywt.idwt(cA, cD, 'db2')
+```
diff --git a/doc/source/regression/modes.rst b/doc/source/regression/modes.rst
deleted file mode 100644
index 9d13ba216..000000000
--- a/doc/source/regression/modes.rst
+++ /dev/null
@@ -1,68 +0,0 @@
-.. _reg-modes:
-
-.. currentmodule:: pywt
-
-
-Signal Extension Modes
-======================
-
-Import :mod:`pywt` first
-
-    >>> import pywt
-
-    >>> def format_array(a):
-    ...     """Consistent array representation across different systems"""
-    ...     import numpy
-    ...     a = numpy.where(numpy.abs(a) < 1e-5, 0, a)
-    ...     return numpy.array2string(a, precision=5, separator=' ', suppress_small=True)
-
-List of available signal extension :ref:`modes <Modes>`:
-
-    >>> print(pywt.Modes.modes)
-    ['zero', 'constant', 'symmetric', 'periodic', 'smooth', 'periodization', 'reflect', 'antisymmetric', 'antireflect']
-
-
-Invalid mode name should rise a :exc:`ValueError`:
-
-    >>> pywt.dwt([1,2,3,4], 'db2', 'invalid')
-    Traceback (most recent call last):
-    ...
-    ValueError: Unknown mode name 'invalid'.
-
-
-You can also refer to modes via :ref:`Modes <Modes>` class attributes:
-
-    >>> x = [1, 2, 1, 5, -1, 8, 4, 6]
-    >>> for mode_name in ['zero', 'constant', 'symmetric', 'reflect', 'periodic', 'smooth', 'periodization']:
-    ...     mode = getattr(pywt.Modes, mode_name)
-    ...     cA, cD = pywt.dwt(x, 'db2', mode)
-    ...     print("Mode: %d (%s)" % (mode, mode_name))
-    Mode: 0 (zero)
-    Mode: 2 (constant)
-    Mode: 1 (symmetric)
-    Mode: 6 (reflect)
-    Mode: 4 (periodic)
-    Mode: 3 (smooth)
-    Mode: 5 (periodization)
-
-
-The default mode is :ref:`symmetric <Modes.symmetric>`:
-
-    >>> cA, cD = pywt.dwt(x, 'db2')
-    >>> print(cA)
-    [ 1.76776695  1.73309178  3.40612438  6.32928585  7.77817459]
-    >>> print(cD)
-    [-0.61237244 -2.15599552 -5.95034847 -1.21545369  1.22474487]
-    >>> print(pywt.idwt(cA, cD, 'db2'))
-    [ 1.  2.  1.  5. -1.  8.  4.  6.]
-
-
-And using a keyword argument:
-
-    >>> cA, cD = pywt.dwt(x, 'db2', mode='symmetric')
-    >>> print(cA)
-    [ 1.76776695  1.73309178  3.40612438  6.32928585  7.77817459]
-    >>> print(cD)
-    [-0.61237244 -2.15599552 -5.95034847 -1.21545369  1.22474487]
-    >>> print(pywt.idwt(cA, cD, 'db2'))
-    [ 1.  2.  1.  5. -1.  8.  4.  6.]
diff --git a/doc/source/regression/multilevel.md b/doc/source/regression/multilevel.md
new file mode 100644
index 000000000..ff0a95dd1
--- /dev/null
+++ b/doc/source/regression/multilevel.md
@@ -0,0 +1,111 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 1.0.0
+    jupytext_version: 1.16.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{include} header.md
+
+```
+
+# Multilevel DWT, IDWT and SWT
+
+## Multilevel DWT decomposition
+
+Here is an example of multilevel DWT decomposition using the `db1` wavelet into
+approximation and detail coefficients:
+
+```{code-cell}
+import pywt
+x = [3, 7, 1, 1, -2, 5, 4, 6]
+db1 = pywt.Wavelet('db1')
+cA3, cD3, cD2, cD1 = pywt.wavedec(x, db1)
+```
+
+```{code-cell}
+cA3
+```
+
+```{code-cell}
+cD3
+```
+
+```{code-cell}
+cD2
+```
+
+```{code-cell}
+cD1
+```
+
+The number of levels in the decomposition can be determined as well:
+
+```{code-cell}
+pywt.dwt_max_level(len(x), db1)
+```
+
+or decompose to a specific level:
+
+```{code-cell}
+cA2, cD2, cD1 = pywt.wavedec(x, db1, mode='constant', level=2)
+```
+
+## Multilevel IDWT reconstruction
+
+```{code-cell}
+coeffs = pywt.wavedec(x, db1)
+pywt.waverec(coeffs, db1)
+```
+
+## Multilevel SWT decomposition
+
+```{code-cell}
+x = [3, 7, 1, 3, -2, 6, 4, 6]
+(cA2, cD2), (cA1, cD1) = pywt.swt(x, db1, level=2)
+```
+
+```{code-cell}
+cA1
+```
+
+```{code-cell}
+cD1
+```
+
+```{code-cell}
+cA2
+```
+
+```{code-cell}
+cD2
+```
+
+```{code-cell}
+[(cA2, cD2)] = pywt.swt(cA1, db1, level=1, start_level=1)
+```
+
+```{code-cell}
+cA2
+```
+
+```{code-cell}
+cD2
+```
+
+```{code-cell}
+coeffs = pywt.swt(x, db1)
+len(coeffs)
+```
+
+```{code-cell}
+pywt.swt_max_level(len(x))
+```
diff --git a/doc/source/regression/multilevel.rst b/doc/source/regression/multilevel.rst
deleted file mode 100644
index c4916834e..000000000
--- a/doc/source/regression/multilevel.rst
+++ /dev/null
@@ -1,64 +0,0 @@
-.. _reg-multilevel:
-
-.. currentmodule:: pywt
-
-Multilevel DWT, IDWT and SWT
-============================
-
-Multilevel DWT decomposition
-----------------------------
-
->>> import pywt
->>> x = [3, 7, 1, 1, -2, 5, 4, 6]
->>> db1 = pywt.Wavelet('db1')
->>> cA3, cD3, cD2, cD1 = pywt.wavedec(x, db1)
->>> print(cA3)
-[ 8.83883476]
->>> print(cD3)
-[-0.35355339]
->>> print(cD2)
-[ 4.  -3.5]
->>> print(cD1)
-[-2.82842712  0.         -4.94974747 -1.41421356]
-
->>> pywt.dwt_max_level(len(x), db1)
-3
-
->>> cA2, cD2, cD1 = pywt.wavedec(x, db1, mode='constant', level=2)
-
-
-Multilevel IDWT reconstruction
-------------------------------
-
->>> coeffs = pywt.wavedec(x, db1)
->>> print(pywt.waverec(coeffs, db1))
-[ 3.  7.  1.  1. -2.  5.  4.  6.]
-
-
-Multilevel SWT decomposition
-----------------------------
-
->>> x = [3, 7, 1, 3, -2, 6, 4, 6]
->>> (cA2, cD2), (cA1, cD1) = pywt.swt(x, db1, level=2)
->>> print(cA1)
-[ 7.07106781  5.65685425  2.82842712  0.70710678  2.82842712  7.07106781
-  7.07106781  6.36396103]
->>> print(cD1)
-[-2.82842712  4.24264069 -1.41421356  3.53553391 -5.65685425  1.41421356
- -1.41421356  2.12132034]
->>> print(cA2)
-[  7.    4.5   4.    5.5   7.    9.5  10.    8.5]
->>> print(cD2)
-[ 3.   3.5  0.  -4.5 -3.   0.5  0.   0.5]
-
->>> [(cA2, cD2)] = pywt.swt(cA1, db1, level=1, start_level=1)
->>> print(cA2)
-[  7.    4.5   4.    5.5   7.    9.5  10.    8.5]
->>> print(cD2)
-[ 3.   3.5  0.  -4.5 -3.   0.5  0.   0.5]
-
->>> coeffs = pywt.swt(x, db1)
->>> len(coeffs)
-3
->>> pywt.swt_max_level(len(x))
-3
diff --git a/doc/source/regression/wavelet.md b/doc/source/regression/wavelet.md
new file mode 100644
index 000000000..688f907c5
--- /dev/null
+++ b/doc/source/regression/wavelet.md
@@ -0,0 +1,235 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 1.0.0
+    jupytext_version: 1.16.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{include} header.md
+
+```
+
+# The Wavelet object
+
+## Wavelet families and builtin Wavelets names
+
+`pywt.Wavelet` objects are really handy carriers of a bunch of DWT-specific
+data like _quadrature mirror filters_ and some general properties associated
+with them.
+
+At first let's go through the methods of creating a `Wavelet` object.
+The easiest and the most convenient way is to use built-in named Wavelets.
+
+These wavelets are organized into groups called wavelet families. The most
+commonly used families are:
+
+```{code-cell}
+import pywt
+pywt.families()
+```
+
+The {func}`wavelist` function with family name passed as an argument is used to
+obtain the list of wavelet names in each family.
+
+```{code-cell}
+for family in pywt.families():
+    print("%s family: " % family + ', '.join(pywt.wavelist(family)))
+```
+
+To get the full list of built-in wavelets' names, just use the `pywt.wavelist` function
+without any arguments.
+
+## Creating Wavelet objects
+
+Now that we know all the names, let's finally create a `Wavelet` object:
+
+```{code-cell}
+w = pywt.Wavelet('db3')
+```
+
+and, that's it!
+
+## Wavelet properties
+
+But what can we do with `Wavelet` objects? Well, they carry some
+interesting pieces of information.
+
+First, let's try printing a `Wavelet` object that we used earlier.
+This shows a brief information about its name, its family name and some
+properties like orthogonality and symmetry.
+
+```{code-cell}
+print(w)
+```
+
+But the most important bits of information are the wavelet filters coefficients,
+which are used in Discrete Wavelet Transform. These coefficients
+can be obtained via the `Wavelet.dec_lo`, `dec_hi`, `rec_lo`,
+and the `rec_hi` attributes, which
+correspond to lowpass & highpass decomposition filters, and lowpass &
+highpass reconstruction filters respectively:
+
+```{code-cell}
+def print_array(arr):
+    print("[%s]" % ", ".join(["%.14f" % x for x in arr]))
+```
+
+Another way to get the filters data is to use the `filter_bank`
+attribute, which returns all four filters in a tuple:
+
+```{code-cell}
+w.filter_bank == (w.dec_lo, w.dec_hi, w.rec_lo, w.rec_hi)
+```
+
+Other properties of a `Wavelet` object are:
+
+`name`, `short_family_name`, and `family_name`:
+
+```{code-cell}
+print(w.name)
+print(w.short_family_name)
+print(w.family_name)
+```
+
+Decomposition (`dec_len`) and reconstruction (`rec_len`) filter lengths:
+
+```{code-cell}
+w.dec_len
+```
+
+```{code-cell}
+w.rec_len
+```
+
+Orthogonality (`orthogonal`) and biorthogonality (`biorthogonal`):
+
+```{code-cell}
+w.orthogonal
+```
+
+```{code-cell}
+w.biorthogonal
+```
+
+Symmetry (`symmetry`):
+
+```{code-cell}
+print(w.symmetry)
+```
+
+Number of vanishing moments for the scaling function `phi` (`vanishing_moments_phi`)
+and the wavelet function `psi` (`vanishing_moments_psi`), associated with the filters:
+
+```{code-cell}
+ w.vanishing_moments_phi
+```
+
+```{code-cell}
+w.vanishing_moments_psi
+```
+
+Now when we know a bit about the builtin Wavelets, let's see how to create
+[custom wavelets](custom-wavelets). These can be created in two ways:
+
+1. Passing the filter bank object that implements the `filter_bank` attribute. The
+   attribute must return four filters coefficients.
+
+```{code-cell}
+class MyHaarFilterBank(object):
+    @property
+    def filter_bank(self):
+        from math import sqrt
+        return (
+          [sqrt(2)/2, sqrt(2)/2], [-sqrt(2)/2, sqrt(2)/2],
+          [sqrt(2)/2, sqrt(2)/2], [sqrt(2)/2, -sqrt(2)/2]
+        )
+
+my_wavelet = pywt.Wavelet('My Haar Wavelet', filter_bank=MyHaarFilterBank())
+```
+
+2. Passing the filters coefficients directly as the `filter_bank` parameter.
+
+```{code-cell}
+from math import sqrt
+my_filter_bank = (
+  [sqrt(2)/2, sqrt(2)/2], [-sqrt(2)/2, sqrt(2)/2],
+  [sqrt(2)/2, sqrt(2)/2], [sqrt(2)/2, -sqrt(2)/2]
+)
+my_wavelet = pywt.Wavelet('My Haar Wavelet', filter_bank=my_filter_bank)
+```
+
+Note that such custom wavelets **will not** have all the properties set
+to correct values and some of them could be missing:
+
+```{code-cell}
+print(my_wavelet)
+```
+
+You can, however, set a couple of them on your own:
+
+```{code-cell}
+my_wavelet.orthogonal = True
+my_wavelet.biorthogonal = True
+```
+
+Let's view the values of the custom wavelet properties again:
+
+```{code-cell}
+print(my_wavelet)
+```
+
+## And now... the `wavefun`!
+
+We all know that the fun with wavelets is in wavelet functions.
+Now, what would this package be without a tool to compute wavelet
+and scaling functions approximations?
+
+This is the purpose of the `wavefun()` method, which is used to
+approximate scaling function (`phi`) and wavelet function (`psi`) at the
+given level of refinement, based on the filters coefficients.
+
+The number of returned values varies depending on the wavelet's
+orthogonality property. For orthogonal wavelets, the result is a tuple
+with the scaling function, wavelet function, and the xgrid coordinates.
+
+```{code-cell}
+w = pywt.Wavelet('sym3')
+w.orthogonal
+```
+
+```{code-cell}
+(phi, psi, x) = w.wavefun(level=5)
+```
+
+For biorthogonal (non-orthogonal) wavelets, different scaling and wavelet
+functions are used for decomposition and reconstruction, and thus, five
+elements are returned: decomposition scaling and wavelet functions
+approximations, reconstruction scaling and wavelet functions approximations,
+and the xgrid.
+
+```{code-cell}
+w = pywt.Wavelet('bior1.3')
+w.orthogonal
+```
+
+```{code-cell}
+(phi_d, psi_d, phi_r, psi_r, x) = w.wavefun(level=5)
+```
+
+:::{seealso}
+You can find live examples of the usage of `wavefun()` and
+images of all the built-in wavelets on the
+[Wavelet Properties Browser](http://wavelets.pybytes.com) page.
+
+However, **this website is no longer actively maintained** and does not
+include every wavelet present in PyWavelets. The precision of the wavelet
+coefficients at that site is also lower than those included in PyWavelets.
+:::
diff --git a/doc/source/regression/wavelet.rst b/doc/source/regression/wavelet.rst
deleted file mode 100644
index e46745371..000000000
--- a/doc/source/regression/wavelet.rst
+++ /dev/null
@@ -1,233 +0,0 @@
-.. _reg-wavelet:
-
-.. currentmodule:: pywt
-
-The Wavelet object
-==================
-
-Wavelet families and builtin Wavelets names
--------------------------------------------
-
-:class:`Wavelet` objects are really a handy carriers of a bunch of DWT-specific
-data like *quadrature mirror filters* and some general properties associated
-with them.
-
-At first let's go through the methods of creating a :class:`Wavelet` object.
-The easiest and the most convenient way is to use builtin named Wavelets.
-
-These wavelets are organized into groups called wavelet families. The most
-commonly used families are:
-
-    >>> import pywt
-    >>> pywt.families()
-    ['haar', 'db', 'sym', 'coif', 'bior', 'rbio', 'dmey', 'gaus', 'mexh', 'morl', 'cgau', 'shan', 'fbsp', 'cmor']
-
-The :func:`wavelist` function with family name passed as an argument is used to
-obtain the list of wavelet names in each family.
-
-    >>> for family in pywt.families():
-    ...     print("%s family: " % family + ', '.join(pywt.wavelist(family)))
-    haar family: haar
-    db family: db1, db2, db3, db4, db5, db6, db7, db8, db9, db10, db11, db12, db13, db14, db15, db16, db17, db18, db19, db20, db21, db22, db23, db24, db25, db26, db27, db28, db29, db30, db31, db32, db33, db34, db35, db36, db37, db38
-    sym family: sym2, sym3, sym4, sym5, sym6, sym7, sym8, sym9, sym10, sym11, sym12, sym13, sym14, sym15, sym16, sym17, sym18, sym19, sym20
-    coif family: coif1, coif2, coif3, coif4, coif5, coif6, coif7, coif8, coif9, coif10, coif11, coif12, coif13, coif14, coif15, coif16, coif17
-    bior family: bior1.1, bior1.3, bior1.5, bior2.2, bior2.4, bior2.6, bior2.8, bior3.1, bior3.3, bior3.5, bior3.7, bior3.9, bior4.4, bior5.5, bior6.8
-    rbio family: rbio1.1, rbio1.3, rbio1.5, rbio2.2, rbio2.4, rbio2.6, rbio2.8, rbio3.1, rbio3.3, rbio3.5, rbio3.7, rbio3.9, rbio4.4, rbio5.5, rbio6.8
-    dmey family: dmey
-    gaus family: gaus1, gaus2, gaus3, gaus4, gaus5, gaus6, gaus7, gaus8
-    mexh family: mexh
-    morl family: morl
-    cgau family: cgau1, cgau2, cgau3, cgau4, cgau5, cgau6, cgau7, cgau8
-    shan family: shan
-    fbsp family: fbsp
-    cmor family: cmor
-
-To get the full list of builtin wavelets' names just use the :func:`wavelist`
-with no argument.
-
-
-Creating Wavelet objects
-------------------------
-
-Now when we know all the names let's finally create a :class:`Wavelet` object:
-
-    >>> w = pywt.Wavelet('db3')
-
-So.. that's it.
-
-
-Wavelet properties
-------------------
-
-But what can we do with :class:`Wavelet` objects? Well, they carry some
-interesting information.
-
-First, let's try printing a :class:`Wavelet` object. This shows a brief
-information about its name, its family name and some properties like
-orthogonality and symmetry.
-
-    >>> print(w)
-    Wavelet db3
-      Family name:    Daubechies
-      Short name:     db
-      Filters length: 6
-      Orthogonal:     True
-      Biorthogonal:   True
-      Symmetry:       asymmetric
-      DWT:            True
-      CWT:            False
-
-But the most important information are the wavelet filters coefficients, which
-are used in :ref:`Discrete Wavelet Transform <ref-dwt>`. These coefficients can
-be obtained via the :attr:`~Wavelet.dec_lo`, :attr:`Wavelet.dec_hi`,
-:attr:`~Wavelet.rec_lo` and :attr:`~Wavelet.rec_hi` attributes, which
-corresponds to lowpass and highpass decomposition filters and lowpass and
-highpass reconstruction filters respectively:
-
-    >>> def print_array(arr):
-    ...     print("[%s]" % ", ".join(["%.14f" % x for x in arr]))
-
-
-
-Another way to get the filters data is to use the :attr:`~Wavelet.filter_bank`
-attribute, which returns all four filters in a tuple:
-
-    >>> w.filter_bank == (w.dec_lo, w.dec_hi, w.rec_lo, w.rec_hi)
-    True
-
-
-Other Wavelet's properties are:
-
-    Wavelet :attr:`~Wavelet.name`, :attr:`~Wavelet.short_family_name` and :attr:`~Wavelet.family_name`:
-
-        >>> print(w.name)
-        db3
-        >>> print(w.short_family_name)
-        db
-        >>> print(w.family_name)
-        Daubechies
-
-    - Decomposition (:attr:`~Wavelet.dec_len`) and reconstruction
-      (:attr:`~.Wavelet.rec_len`) filter lengths:
-
-        >>> int(w.dec_len) # int() is for normalizing longs and ints for doctest
-        6
-        >>> int(w.rec_len)
-        6
-
-    - Orthogonality (:attr:`~Wavelet.orthogonal`) and biorthogonality (:attr:`~Wavelet.biorthogonal`):
-
-        >>> w.orthogonal
-        True
-        >>> w.biorthogonal
-        True
-
-    - Symmetry (:attr:`~Wavelet.symmetry`):
-
-        >>> print(w.symmetry)
-        asymmetric
-
-    - Number of vanishing moments for the scaling function ``phi``
-      (:attr:`~Wavelet.vanishing_moments_phi`) and the wavelet function ``psi``
-      (:attr:`~Wavelet.vanishing_moments_psi`) associated with the filters:
-
-        >>> w.vanishing_moments_phi
-        0
-        >>> w.vanishing_moments_psi
-        3
-
-Now when we know a bit about the builtin Wavelets, let's see how to create
-:ref:`custom Wavelets <custom-wavelets>` objects. These can be done in two ways:
-
-    1) Passing the filter bank object that implements the ``filter_bank``
-       attribute. The attribute must return four filters coefficients.
-
-       >>> class MyHaarFilterBank(object):
-       ...     @property
-       ...     def filter_bank(self):
-       ...         from math import sqrt
-       ...         return ([sqrt(2)/2, sqrt(2)/2], [-sqrt(2)/2, sqrt(2)/2],
-       ...                 [sqrt(2)/2, sqrt(2)/2], [sqrt(2)/2, -sqrt(2)/2])
-
-
-       >>> my_wavelet = pywt.Wavelet('My Haar Wavelet', filter_bank=MyHaarFilterBank())
-
-
-    2) Passing the filters coefficients directly as the ``filter_bank``
-       parameter.
-
-        >>> from math import sqrt
-        >>> my_filter_bank = ([sqrt(2)/2, sqrt(2)/2], [-sqrt(2)/2, sqrt(2)/2],
-        ...                   [sqrt(2)/2, sqrt(2)/2], [sqrt(2)/2, -sqrt(2)/2])
-        >>> my_wavelet = pywt.Wavelet('My Haar Wavelet', filter_bank=my_filter_bank)
-
-
-Note that such custom wavelets **will not** have all the properties set
-to correct values:
-
-    >>> print(my_wavelet)
-    Wavelet My Haar Wavelet
-      Family name:
-      Short name:
-      Filters length: 2
-      Orthogonal:     False
-      Biorthogonal:   False
-      Symmetry:       unknown
-      DWT:            True
-      CWT:            False
-
-    You can however set a couple of them on your own:
-
-    >>> my_wavelet.orthogonal = True
-    >>> my_wavelet.biorthogonal = True
-
-    >>> print(my_wavelet)
-    Wavelet My Haar Wavelet
-      Family name:
-      Short name:
-      Filters length: 2
-      Orthogonal:     True
-      Biorthogonal:   True
-      Symmetry:       unknown
-      DWT:            True
-      CWT:            False
-
-And now... the ``wavefun``!
----------------------------
-
-We all know that the fun with wavelets is in wavelet functions.
-Now what would be this package without a tool to compute wavelet
-and scaling functions approximations?
-
-This is the purpose of the :meth:`~Wavelet.wavefun` method, which is used to
-approximate scaling function (``phi``) and wavelet function (``psi``) at the
-given level of refinement, based on the filters coefficients.
-
-The number of returned values varies depending on the wavelet's
-orthogonality property. For orthogonal wavelets the result is tuple
-with scaling function, wavelet function and xgrid coordinates.
-
-    >>> w = pywt.Wavelet('sym3')
-    >>> w.orthogonal
-    True
-    >>> (phi, psi, x) = w.wavefun(level=5)
-
-For biorthogonal (non-orthogonal) wavelets different scaling and wavelet
-functions are used for decomposition and reconstruction, and thus five
-elements are returned: decomposition scaling and wavelet functions
-approximations, reconstruction scaling and wavelet functions approximations,
-and the xgrid.
-
-    >>> w = pywt.Wavelet('bior1.3')
-    >>> w.orthogonal
-    False
-    >>> (phi_d, psi_d, phi_r, psi_r, x) = w.wavefun(level=5)
-
-.. seealso::
-      You can find live examples of :meth:`~Wavelet.wavefun` usage and
-      images of all the built-in wavelets on the
-      `Wavelet Properties Browser <http://wavelets.pybytes.com>`_ page.
-      However, **this website is no longer actively maintained** and does not
-      include every wavelet present in PyWavelets. The precision of the wavelet
-      coefficients at that site is also lower than those included in
-      PyWavelets.
diff --git a/doc/source/regression/wp.md b/doc/source/regression/wp.md
new file mode 100644
index 000000000..b3200a80a
--- /dev/null
+++ b/doc/source/regression/wp.md
@@ -0,0 +1,382 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 1.0.0
+    jupytext_version: 1.16.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{include} header.md
+
+```
+
+# Wavelet Packets
+
+## `import pywt` and construct a helper function
+
+```{code-cell}
+import pywt
+```
+
+This helper function can format arrays in a consistent manner across
+different systems. Please note that this function is just for the purpose of
+this example and is not part of the PyWavelets library, and it is not necessary
+or required to use it in your own code:
+
+```{code-cell}
+def format_array(a):
+    """Consistent array representation across different systems"""
+    import numpy
+    a = numpy.where(numpy.abs(a) < 1e-5, 0, a)
+    return numpy.array2string(a, precision=5, separator=' ', suppress_small=True)
+```
+
+## Create Wavelet Packet structure
+
+Okay, let's create a sample instance of the `WaveletPacket` class:
+
+```{code-cell}
+x = [1, 2, 3, 4, 5, 6, 7, 8]
+wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
+```
+
+The input data and decomposition coefficients are stored in the
+`WaveletPacket.data` attribute:
+
+```{code-cell}
+print(wp.data)
+```
+
+Nodes are identified by their paths, i.e., `Node.path`. For the root
+node, the path is `''` and the decomposition level is `0`.
+
+```{code-cell}
+print(repr(wp.path))
+```
+
+```{code-cell}
+print(wp.level)
+```
+
+The `maxlevel` attribute, if not given as param in the constructor, is
+automatically computed:
+
+```{code-cell}
+print(wp['ad'].maxlevel)
+```
+
+## Traversing the WP tree
+
+### Accessing subnodes
+
+```{code-cell}
+x = [1, 2, 3, 4, 5, 6, 7, 8]
+wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
+```
+
+First check what is the maximum level of decomposition:
+
+```{code-cell}
+print(wp.maxlevel)
+```
+
+and try accessing subnodes of the WP tree.
+
+1st level:
+
+```{code-cell}
+print(wp['a'].data)
+```
+
+```{code-cell}
+print(wp['a'].path)
+```
+
+2nd level:
+
+```{code-cell}
+print(wp['aa'].data)
+```
+
+```{code-cell}
+print(wp['aa'].path)
+```
+
+3rd level:
+
+```{code-cell}
+print(wp['aaa'].data)
+```
+
+```{code-cell}
+print(wp['aaa'].path)
+```
+
+Oops, we have reached the maximum level of decomposition and received an `IndexError`:
+
+```{code-cell}
+---
+tags: [raises-exception, remove-output]
+---
+print(wp['aaaa'].data)
+```
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{code-block} python
+:class: pywt-handcoded-cell-output
+Traceback (most recent call last):
+...
+IndexError: Path length is out of range.
+```
+
+Now, try an invalid path:
+
+```{code-cell}
+---
+tags: [raises-exception, remove-output]
+---
+print(wp['ac'])
+```
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{code-block} python
+:class: pywt-handcoded-cell-output
+Traceback (most recent call last):
+...
+ValueError: Subnode name must be in ['a', 'd'], not 'c'.
+```
+
+which just yielded a `ValueError`.
+
+### Accessing `Node`'s attributes
+
+The `WaveletPacket` object is a tree data structure, which evaluates to a set
+of `Node` objects. `WaveletPacket` is just a special subclass
+of the `Node` class (which in turn inherits from the `BaseNode` class).
+
+Tree nodes can be accessed using the `obj[x]` (`Node.__getitem__`)
+operator.
+
+Each tree node has a set of attributes: `data`, `path`,
+`node_name`, `parent`, `level`,
+`maxlevel` and `mode`.
+
+```{code-cell}
+x = [1, 2, 3, 4, 5, 6, 7, 8]
+wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
+```
+
+```{code-cell}
+print(wp['ad'].data)
+```
+
+```{code-cell}
+print(wp['ad'].path)
+```
+
+```{code-cell}
+print(wp['ad'].node_name)
+```
+
+```{code-cell}
+print(wp['ad'].parent.path)
+```
+
+```{code-cell}
+print(wp['ad'].level)
+```
+
+```{code-cell}
+print(wp['ad'].maxlevel)
+```
+
+```{code-cell}
+print(wp['ad'].mode)
+```
+
+### Collecting nodes
+
+```{code-cell}
+x = [1, 2, 3, 4, 5, 6, 7, 8]
+wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
+```
+
+We can get all nodes on the particular level either in `natural` order:
+
+```{code-cell}
+print([node.path for node in wp.get_level(3, 'natural')])
+```
+
+or sorted based on the band frequency (`freq`):
+
+```{code-cell}
+print([node.path for node in wp.get_level(3, 'freq')])
+```
+
+Note that `WaveletPacket.get_level()` also performs automatic decomposition
+until it reaches the specified `level`.
+
+## Reconstructing data from Wavelet Packets
+
+```{code-cell}
+x = [1, 2, 3, 4, 5, 6, 7, 8]
+wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
+```
+
+Now create a new instance of the `WaveletPacket` class and set its nodes
+with some data.
+
+```{code-cell}
+new_wp = pywt.WaveletPacket(data=None, wavelet='db1', mode='symmetric')
+```
+
+```{code-cell}
+new_wp['aa'] = wp['aa'].data
+new_wp['ad'] = [-2., -2.]
+```
+
+For convenience, `Node.data` gets automatically extracted from the `Node` object:
+
+```{code-cell}
+new_wp['d'] = wp['d']
+```
+
+And reconstruct the data from the `aa`, `ad` and `d` packets.
+
+```{code-cell}
+print(new_wp.reconstruct(update=False))
+```
+
+If the `update` param in the reconstruct method is set to `False`, the
+node's `Node.data` will not be updated.
+
+```{code-cell}
+print(new_wp.data)
+```
+
+Otherwise, the `Node.data` attribute will be set to the reconstructed
+value.
+
+```{code-cell}
+print(new_wp.reconstruct(update=True))
+```
+
+```{code-cell}
+print(new_wp.data)
+```
+
+```{code-cell}
+print([n.path for n in new_wp.get_leaf_nodes(False)])
+```
+
+```{code-cell}
+print([n.path for n in new_wp.get_leaf_nodes(True)])
+```
+
+## Removing nodes from Wavelet Packet tree
+
+Let's create some sample data:
+
+```{code-cell}
+x = [1, 2, 3, 4, 5, 6, 7, 8]
+wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
+```
+
+First, let's start with a tree decomposition at level 2. The leaf nodes in the tree are:
+
+```{code-cell}
+dummy = wp.get_level(2)
+for n in wp.get_leaf_nodes(False):
+    print(n.path, format_array(n.data))
+```
+
+```{code-cell}
+node = wp['ad']
+print(node)
+```
+
+To remove a node from the WP tree, use Python's `del obj[x]` (`Node.__delitem__`):
+
+```{code-cell}
+del wp['ad']
+```
+
+The leaf nodes that left in the tree are:
+
+```{code-cell}
+for n in wp.get_leaf_nodes():
+    print(n.path, format_array(n.data))
+```
+
+And the reconstruction is:
+
+```{code-cell}
+print(wp.reconstruct())
+```
+
+Now, restore the deleted node value.
+
+```{code-cell}
+wp['ad'].data = node.data
+```
+
+Printing leaf nodes and tree reconstruction confirms the original state of the
+tree:
+
+```{code-cell}
+for n in wp.get_leaf_nodes(False):
+    print(n.path, format_array(n.data))
+```
+
+```{code-cell}
+print(wp.reconstruct())
+```
+
+## Lazy evaluation
+
+:::{note}
+This section is for demonstration of PyWavelets' internal purposes
+only. Do not rely on the attribute access to nodes as presented in
+this example.
+:::
+
+```{code-cell}
+x = [1, 2, 3, 4, 5, 6, 7, 8]
+wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
+```
+
+At first the wp's attribute `a` is `None`:
+
+```{code-cell}
+print(wp.a)
+```
+
+**Remember that you should not rely on the attribute access.**
+
+At the first attempt to access the node, it is computed via the decomposition
+of its parent node (which is the `wp` object itself).
+
+```{code-cell}
+print(wp['a'])
+```
+
+Now `wp.a` is set to the newly created node:
+
+```{code-cell}
+print(wp.a)
+```
+
+And so is `wp.d`:
+
+```{code-cell}
+print(wp.d)
+```
diff --git a/doc/source/regression/wp.rst b/doc/source/regression/wp.rst
deleted file mode 100644
index 8afc12ec4..000000000
--- a/doc/source/regression/wp.rst
+++ /dev/null
@@ -1,303 +0,0 @@
-.. _reg-wp:
-
-.. currentmodule:: pywt
-
->>> from __future__ import print_function
-
-Wavelet Packets
-===============
-
-Import pywt
------------
-
-    >>> import pywt
-
-    >>> def format_array(a):
-    ...     """Consistent array representation across different systems"""
-    ...     import numpy
-    ...     a = numpy.where(numpy.abs(a) < 1e-5, 0, a)
-    ...     return numpy.array2string(a, precision=5, separator=' ', suppress_small=True)
-
-
-Create Wavelet Packet structure
--------------------------------
-
-Ok, let's create a sample :class:`WaveletPacket`:
-
-    >>> x = [1, 2, 3, 4, 5, 6, 7, 8]
-    >>> wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
-
-The input ``data`` and decomposition coefficients are stored in the
-:attr:`WaveletPacket.data` attribute:
-
-    >>> print(wp.data)
-    [1, 2, 3, 4, 5, 6, 7, 8]
-
-:class:`Nodes <Node>` are identified by :attr:`paths <~Node.path>`. For the root
-node the path is ``''`` and the decomposition level is ``0``.
-
-    >>> print(repr(wp.path))
-    ''
-    >>> print(wp.level)
-    0
-
-The ``maxlevel``, if not given as param in the constructor, is automatically
-computed:
-
-    >>> print(wp['ad'].maxlevel)
-    3
-
-
-Traversing WP tree:
--------------------
-
-Accessing subnodes:
-~~~~~~~~~~~~~~~~~~~
-
->>> x = [1, 2, 3, 4, 5, 6, 7, 8]
->>> wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
-
-First check what is the maximum level of decomposition:
-
-    >>> print(wp.maxlevel)
-    3
-
-and try accessing subnodes of the WP tree:
-
-    * 1st level:
-
-        >>> print(wp['a'].data)
-        [  2.12132034   4.94974747   7.77817459  10.60660172]
-        >>> print(wp['a'].path)
-        a
-
-    * 2nd level:
-
-        >>> print(wp['aa'].data)
-        [  5.  13.]
-        >>> print(wp['aa'].path)
-        aa
-
-
-    * 3rd level:
-
-        >>> print(wp['aaa'].data)
-        [ 12.72792206]
-        >>> print(wp['aaa'].path)
-        aaa
-
-
-      Ups, we have reached the maximum level of decomposition and got an
-      :exc:`IndexError`:
-
-        >>> print(wp['aaaa'].data)
-        Traceback (most recent call last):
-        ...
-        IndexError: Path length is out of range.
-
-Now try some invalid path:
-
-    >>> print(wp['ac'])
-    Traceback (most recent call last):
-    ...
-    ValueError: Subnode name must be in ['a', 'd'], not 'c'.
-
-which just yielded a :exc:`ValueError`.
-
-
-Accessing Node's attributes:
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-:class:`WaveletPacket` object is a tree data structure, which evaluates to a set
-of :class:`Node` objects. :class:`WaveletPacket` is just a special subclass
-of the :class:`Node` class (which in turn inherits from the :class:`BaseNode`).
-
-Tree nodes can be accessed using the ``obj[x]`` (:meth:`Node.__getitem__`)
-operator.
-Each tree node has a set of attributes: :attr:`~Node.data`, :attr:`~Node.path`,
-:attr:`~Node.node_name`, :attr:`~Node.parent`, :attr:`~Node.level`,
-:attr:`~Node.maxlevel` and :attr:`~Node.mode`.
-
->>> x = [1, 2, 3, 4, 5, 6, 7, 8]
->>> wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
-
->>> print(wp['ad'].data)
-[-2. -2.]
-
->>> print(wp['ad'].path)
-ad
-
->>> print(wp['ad'].node_name)
-d
-
->>> print(wp['ad'].parent.path)
-a
-
->>> print(wp['ad'].level)
-2
-
->>> print(wp['ad'].maxlevel)
-3
-
->>> print(wp['ad'].mode)
-symmetric
-
-
-Collecting nodes
-~~~~~~~~~~~~~~~~
-
->>> x = [1, 2, 3, 4, 5, 6, 7, 8]
->>> wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
-
-
-We can get all nodes on the particular level either in ``natural`` order:
-
-    >>> print([node.path for node in wp.get_level(3, 'natural')])
-    ['aaa', 'aad', 'ada', 'add', 'daa', 'dad', 'dda', 'ddd']
-
-or sorted based on the band frequency (``freq``):
-
-    >>> print([node.path for node in wp.get_level(3, 'freq')])
-    ['aaa', 'aad', 'add', 'ada', 'dda', 'ddd', 'dad', 'daa']
-
-Note that :meth:`WaveletPacket.get_level` also performs automatic decomposition
-until it reaches the specified ``level``.
-
-
-Reconstructing data from Wavelet Packets:
------------------------------------------
-
->>> x = [1, 2, 3, 4, 5, 6, 7, 8]
->>> wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
-
-
-Now create a new :class:`Wavelet Packet <WaveletPacket>` and set its nodes with
-some data.
-
-    >>> new_wp = pywt.WaveletPacket(data=None, wavelet='db1', mode='symmetric')
-
-    >>> new_wp['aa'] = wp['aa'].data
-    >>> new_wp['ad'] = [-2., -2.]
-
-For convenience, :attr:`Node.data` gets automatically extracted from the
-:class:`Node` object:
-
-    >>> new_wp['d'] = wp['d']
-
-And reconstruct the data from the ``aa``, ``ad`` and ``d`` packets.
-
-    >>> print(new_wp.reconstruct(update=False))
-    [ 1.  2.  3.  4.  5.  6.  7.  8.]
-
-If the ``update`` param in the reconstruct method is set to ``False``, the
-node's :attr:`~Node.data` will not be updated.
-
-    >>> print(new_wp.data)
-    None
-
-Otherwise, the :attr:`~Node.data` attribute will be set to the reconstructed
-value.
-
-    >>> print(new_wp.reconstruct(update=True))
-    [ 1.  2.  3.  4.  5.  6.  7.  8.]
-    >>> print(new_wp.data)
-    [ 1.  2.  3.  4.  5.  6.  7.  8.]
-
-
->>> print([n.path for n in new_wp.get_leaf_nodes(False)])
-['aa', 'ad', 'd']
-
->>> print([n.path for n in new_wp.get_leaf_nodes(True)])
-['aaa', 'aad', 'ada', 'add', 'daa', 'dad', 'dda', 'ddd']
-
-
-Removing nodes from Wavelet Packet tree:
-----------------------------------------
-
-Let's create a sample data:
-
-    >>> x = [1, 2, 3, 4, 5, 6, 7, 8]
-    >>> wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
-
-First, start with a tree decomposition at level 2. Leaf nodes in the tree are:
-
-    >>> dummy = wp.get_level(2)
-    >>> for n in wp.get_leaf_nodes(False):
-    ...     print(n.path, format_array(n.data))
-    aa [  5.  13.]
-    ad [-2. -2.]
-    da [-1. -1.]
-    dd [ 0.  0.]
-
-    >>> node = wp['ad']
-    >>> print(node)
-    ad: [-2. -2.]
-
-To remove a node from the WP tree, use Python's ``del obj[x]``
-(:class:`Node.__delitem__`):
-
-    >>> del wp['ad']
-
-The leaf nodes that left in the tree are:
-
-    >>> for n in wp.get_leaf_nodes():
-    ...     print(n.path, format_array(n.data))
-    aa [  5.  13.]
-    da [-1. -1.]
-    dd [ 0.  0.]
-
-And the reconstruction is:
-
-    >>> print(wp.reconstruct())
-    [ 2.  3.  2.  3.  6.  7.  6.  7.]
-
-Now restore the deleted node value.
-
-    >>> wp['ad'].data = node.data
-
-Printing leaf nodes and tree reconstruction confirms the original state of the
-tree:
-
-    >>> for n in wp.get_leaf_nodes(False):
-    ...     print(n.path, format_array(n.data))
-    aa [  5.  13.]
-    ad [-2. -2.]
-    da [-1. -1.]
-    dd [ 0.  0.]
-
-    >>> print(wp.reconstruct())
-    [ 1.  2.  3.  4.  5.  6.  7.  8.]
-
-
-Lazy evaluation:
-----------------
-
-.. note:: This section is for demonstration of pywt internals purposes
-          only. Do not rely on the attribute access to nodes as presented in
-          this example.
-
->>> x = [1, 2, 3, 4, 5, 6, 7, 8]
->>> wp = pywt.WaveletPacket(data=x, wavelet='db1', mode='symmetric')
-
-1) At first the wp's attribute ``a`` is None
-
-   >>> print(wp.a)
-   None
-
-   **Remember that you should not rely on the attribute access.**
-
-2) At first attempt to access the node it is computed via decomposition
-   of its parent node (the wp object itself).
-
-   >>> print(wp['a'])
-   a: [  2.12132034   4.94974747   7.77817459  10.60660172]
-
-3) Now the ``wp.a`` is set to the newly created node:
-
-   >>> print(wp.a)
-   a: [  2.12132034   4.94974747   7.77817459  10.60660172]
-
-   And so is ``wp.d``:
-
-   >>> print(wp.d)
-   d: [-0.70710678 -0.70710678 -0.70710678 -0.70710678]
diff --git a/doc/source/regression/wp2d.md b/doc/source/regression/wp2d.md
new file mode 100644
index 000000000..178d1e689
--- /dev/null
+++ b/doc/source/regression/wp2d.md
@@ -0,0 +1,430 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 1.0.0
+    jupytext_version: 1.16.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{include} header.md
+
+```
+
+# 2D Wavelet Packets
+
+## Import pywt
+
+```{code-cell}
+import pywt
+import numpy
+```
+
+## Create 2D Wavelet Packet structure
+
+Start with preparing test data:
+
+```{code-cell}
+x = numpy.array([[1, 2, 3, 4, 5, 6, 7, 8]] * 8, 'd')
+print(x)
+```
+
+Now create a 2D [Wavelet Packet](ref-wp) object:
+
+```{code-cell}
+wp = pywt.WaveletPacket2D(data=x, wavelet='db1', mode='symmetric')
+```
+
+The input `data` and decomposition coefficients are stored in the
+`WaveletPacket2D.data` attribute:
+
+```{code-cell}
+print(wp.data)
+```
+
+Nodes (the `Node2D` class) are identified by paths. For the root node, the path is
+`''` and the decomposition level is `0`.
+
+```{code-cell}
+print(repr(wp.path))
+```
+
+```{code-cell}
+print(wp.level)
+```
+
+`WaveletPacket2D.maxlevel`, if not given in the constructor, is
+automatically computed based on the data size:
+
+```{code-cell}
+print(wp.maxlevel)
+```
+
+## Traversing WP tree
+
+Wavelet Packet nodes (`Node2D`) are arranged in a tree. Each node in a WP
+tree is uniquely identified and addressed by a `Node2D.path` string.
+
+In the 1D `WaveletPacket` case nodes were accessed using `'a'`
+(approximation) and `'d'` (details) path names (each node has two 1D
+children).
+
+Because now we deal with a bit more complex structure (each node has four
+children), we have four basic path names based on the dwt 2D output convention
+to address the WP2D structure:
+
+- `a` - LL, low-low coefficients
+- `h` - LH, low-high coefficients
+- `v` - HL, high-low coefficients
+- `d` - HH, high-high coefficients
+
+In other words, subnode naming corresponds to the `dwt2` function output
+naming convention (as wavelet packet transform is based on the dwt2 transform):
+
+```
+                            -------------------
+                            |        |        |
+                            | cA(LL) | cH(LH) |
+                            |        |        |
+(cA, (cH, cV, cD))  <--->   -------------------
+                            |        |        |
+                            | cV(HL) | cD(HH) |
+                            |        |        |
+                            -------------------
+
+   (fig.1: DWT 2D output and interpretation)
+```
+
+Knowing what the nodes names are, we can now access them using the indexing
+operator `obj[x]` (`WaveletPacket2D.__getitem__`):
+
+```{code-cell}
+print(wp['a'].data)
+```
+
+```{code-cell}
+print(wp['h'].data)
+```
+
+```{code-cell}
+print(wp['v'].data)
+```
+
+```{code-cell}
+print(wp['d'].data)
+```
+
+Similarly, a subnode of a subnode can be accessed by:
+
+```{code-cell}
+print(wp['aa'].data)
+```
+
+Indexing base 2D (`WaveletPacket2D`) (as well as 1D `WaveletPacket`)
+using compound paths is just the same as indexing the WP subnode:
+
+```{code-cell}
+node = wp['a']
+print(node['a'].data)
+```
+
+```{code-cell}
+print(wp['a']['a'].data is wp['aa'].data)
+```
+
+Following down the decomposition path:
+
+```{code-cell}
+print(wp['aaa'].data)
+```
+
+```{code-cell}
+---
+tags: [raises-exception, remove-output]
+---
+print(wp['aaaa'].data)
+```
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{code-block} python
+:class: pywt-handcoded-cell-output
+Traceback (most recent call last):
+...
+IndexError: Path length is out of range.
+```
+
+Oops, we have reached the maximum level of decomposition for the `'aaaa'` path,
+which, by the way, was:
+
+```{code-cell}
+print(wp.maxlevel)
+```
+
+Now, try an invalid path:
+
+```{code-cell}
+---
+tags: [raises-exception, remove-output]
+---
+print(wp['f'])
+```
+
++++ {"tags": ["jupyterlite_sphinx_strip"]}
+
+```{code-block} python
+:class: pywt-handcoded-cell-output
+Traceback (most recent call last):
+...
+ValueError: Subnode name must be in ['a', 'h', 'v', 'd'], not 'f'.
+```
+
+### Accessing Node2D's attributes
+
+`WaveletPacket2D` is a tree data structure, which evaluates to a set
+of `Node2D` objects. `WaveletPacket2D` is just a special the `Node2D` class (which in turn inherits from a `BaseNode` class
+just like with `Node` and `WaveletPacket` for the 1D case).
+
+```{code-cell}
+print(wp['av'].data)
+```
+
+```{code-cell}
+print(wp['av'].path)
+```
+
+```{code-cell}
+print(wp['av'].node_name)
+```
+
+```{code-cell}
+print(wp['av'].parent.path)
+```
+
+```{code-cell}
+print(wp['av'].parent.data)
+```
+
+```{code-cell}
+print(wp['av'].level)
+```
+
+```{code-cell}
+print(wp['av'].maxlevel)
+```
+
+```{code-cell}
+print(wp['av'].mode)
+```
+
+### Collecting nodes
+
+We can get all nodes on the particular level using the
+`WaveletPacket2D.get_level` method:
+
+0 level - the root `wp` node:
+
+```{code-cell}
+len(wp.get_level(0))
+```
+
+```{code-cell}
+print([node.path for node in wp.get_level(0)])
+```
+
+- 1st level of decomposition:
+
+```{code-cell}
+len(wp.get_level(1))
+```
+
+```{code-cell}
+print([node.path for node in wp.get_level(1)])
+```
+
+2nd level of decomposition:
+
+```{code-cell}
+len(wp.get_level(2))
+```
+
+```{code-cell}
+paths = [node.path for node in wp.get_level(2)]
+for i, path in enumerate(paths):
+    if (i+1) % 4 == 0:
+        print(path)
+    else:
+        print(path, end=' ')
+```
+
+3rd level of decomposition:
+
+```{code-cell}
+print(len(wp.get_level(3)))
+```
+
+```{code-cell}
+paths = [node.path for node in wp.get_level(3)]
+for i, path in enumerate(paths):
+    if (i+1) % 8 == 0:
+        print(path)
+    else:
+        print(path, end=' ')
+```
+
+Note that `WaveletPacket2D.get_level` performs automatic decomposition
+until it reaches the given level.
+
+## Reconstructing data from Wavelet Packets
+
+Let's create a new empty 2D Wavelet Packet structure and set its nodes
+values with known data from the previous examples:
+
+```{code-cell}
+new_wp = pywt.WaveletPacket2D(data=None, wavelet='db1', mode='symmetric')
+```
+
+```{code-cell}
+new_wp['vh'] = wp['vh'].data  # [[0.0, 0.0], [0.0, 0.0]]
+new_wp['vv'] = wp['vh'].data  # [[0.0, 0.0], [0.0, 0.0]]
+new_wp['vd'] = [[0.0, 0.0], [0.0, 0.0]]
+```
+
+```{code-cell}
+new_wp['a'] = [[3.0, 7.0, 11.0, 15.0], [3.0, 7.0, 11.0, 15.0],
+              [3.0, 7.0, 11.0, 15.0], [3.0, 7.0, 11.0, 15.0]]
+new_wp['d'] = [[0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0],
+              [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0]]
+```
+
+For convenience, `Node2D.data` gets automatically extracted from the
+base `Node2D` object:
+
+```{code-cell}
+new_wp['h'] = wp['h'] # all zeros
+```
+
+Note: just remember to not assign to the `node.data parameter directly (TODO).
+
+And reconstruct the data from the `a`, `d`, `vh`, `vv`, `vd` and `h`
+packets (Note that `va` node was not set and the WP tree is "not complete"
+\- the `va` branch will be treated as _zero-array_):
+
+```{code-cell}
+print(new_wp.reconstruct(update=False))
+```
+
+Now set the `va` node with the known values and do the reconstruction again:
+
+```{code-cell}
+new_wp['va'] = wp['va'].data # [[-2.0, -2.0], [-2.0, -2.0]]
+print(new_wp.reconstruct(update=False))
+```
+
+which is just the same as the base sample data `x`.
+
+Of course we can go the other way and remove nodes from the tree. If we delete
+the `va` node, again, we get the "not complete" tree from one of the previous
+examples:
+
+```{code-cell}
+del new_wp['va']
+print(new_wp.reconstruct(update=False))
+```
+
+Just restore the node before the next examples:
+
+```{code-cell}
+new_wp['va'] = wp['va'].data
+```
+
+If the `update` param in the `WaveletPacket2D.reconstruct` method is set
+to `False`, the node's `Node2D.data` attribute will not be updated.
+
+```{code-cell}
+print(new_wp.data)
+```
+
+Otherwise, the `WaveletPacket2D.data` attribute will be set to the
+reconstructed value.
+
+```{code-cell}
+print(new_wp.reconstruct(update=True))
+```
+
+```{code-cell}
+print(new_wp.data)
+```
+
+Since we have an interesting WP structure built, it is a good occasion to
+present the `WaveletPacket2D.get_leaf_nodes()` method, which collects
+non-zero leaf nodes from the WP tree:
+
+```{code-cell}
+print([n.path for n in new_wp.get_leaf_nodes()])
+```
+
+Passing the `decompose = True` parameter to the method will force the WP
+object to do a full decomposition up to the _maximum level_ of decomposition:
+
+```{code-cell}
+paths = [n.path for n in new_wp.get_leaf_nodes(decompose=True)]
+len(paths)
+```
+
+```{code-cell}
+for i, path in enumerate(paths):
+    if (i+1) % 8 == 0:
+        print(path)
+    else:
+        try:
+            print(path, end=' ')
+        except:
+            print(path, end=' ')
+```
+
+## Lazy evaluation
+
+:::{note}
+This section is for the demonstration of PyWavelets' internals' purposes
+only. Do not rely on the attribute access to nodes as presented in
+this example.
+:::
+
+```{code-cell}
+x = numpy.array([[1, 2, 3, 4, 5, 6, 7, 8]] * 8)
+wp = pywt.WaveletPacket2D(data=x, wavelet='db1', mode='symmetric')
+```
+
+At first, the `wp`'s attribute `a` is `None`
+
+```{code-cell}
+print(wp.a)
+```
+
+**Remember that you should not rely on the attribute access.**
+
+During the first attempt to access the node it is computed
+via decomposition of its parent node (the wp object itself).
+
+```{code-cell}
+print(wp['a'])
+```
+
+Now the `a` is set to the newly created node:
+
+```{code-cell}
+print(wp.a)
+```
+
+And so is `wp.d`:
+
+```{code-cell}
+print(wp.d)
+```
diff --git a/doc/source/regression/wp2d.rst b/doc/source/regression/wp2d.rst
deleted file mode 100644
index 8d70d9185..000000000
--- a/doc/source/regression/wp2d.rst
+++ /dev/null
@@ -1,429 +0,0 @@
-.. _reg-wp2d:
-
-.. currentmodule:: pywt
-
-
-2D Wavelet Packets
-==================
-
-Import pywt
------------
-
->>> from __future__ import print_function
->>> import pywt
->>> import numpy
-
-
-Create 2D Wavelet Packet structure
-----------------------------------
-
-Start with preparing test data:
-
-    >>> x = numpy.array([[1, 2, 3, 4, 5, 6, 7, 8]] * 8, 'd')
-    >>> print(x)
-    [[ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]]
-
-Now create a :class:`2D Wavelet Packet <WaveletPacket2D>` object:
-
-    >>> wp = pywt.WaveletPacket2D(data=x, wavelet='db1', mode='symmetric')
-
-The input ``data`` and decomposition coefficients are stored in the
-:attr:`WaveletPacket2D.data` attribute:
-
-    >>> print(wp.data)
-    [[ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]]
-
-:class:`Nodes <Node2D>` are identified by paths. For the root node the path is
-``''`` and the decomposition level is ``0``.
-
-    >>> print(repr(wp.path))
-    ''
-    >>> print(wp.level)
-    0
-
-The :attr:`WaveletPacket2D.maxlevel`, if not given in the constructor, is
-automatically computed based on the data size:
-
-    >>> print(wp.maxlevel)
-    3
-
-
-Traversing WP tree:
--------------------
-
-Wavelet Packet :class:`nodes <Node2D>` are arranged in a tree. Each node in a WP
-tree is uniquely identified and addressed by a :attr:`~Node2D.path` string.
-
-In the 1D :class:`WaveletPacket` case nodes were accessed using ``'a'``
-(approximation) and ``'d'`` (details) path names (each node has two 1D
-children).
-
-Because now we deal with a bit more complex structure (each node has four
-children), we have four basic path names based on the dwt 2D output convention
-to address the WP2D structure:
-
-    * ``a`` - LL, low-low coefficients
-    * ``h`` - LH, low-high coefficients
-    * ``v`` - HL, high-low coefficients
-    * ``d`` - HH, high-high coefficients
-
-In other words, subnode naming corresponds to the :func:`dwt2` function output
-naming convention (as wavelet packet transform is based on the dwt2 transform)::
-
-                                -------------------
-                                |        |        |
-                                | cA(LL) | cH(LH) |
-                                |        |        |
-    (cA, (cH, cV, cD))  <--->   -------------------
-                                |        |        |
-                                | cV(HL) | cD(HH) |
-                                |        |        |
-                                -------------------
-
-       (fig.1: DWT 2D output and interpretation)
-
-
-Knowing what the nodes names are, we can now access them using the indexing
-operator ``obj[x]`` (:meth:`WaveletPacket2D.__getitem__`):
-
-    >>> print(wp['a'].data)
-    [[  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]]
-    >>> print(wp['h'].data)
-    [[ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]]
-    >>> print(wp['v'].data)
-    [[-1. -1. -1. -1.]
-     [-1. -1. -1. -1.]
-     [-1. -1. -1. -1.]
-     [-1. -1. -1. -1.]]
-    >>> print(wp['d'].data)
-    [[ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]]
-
-Similarly, a subnode of a subnode can be accessed by:
-
-    >>> print(wp['aa'].data)
-    [[ 10.  26.]
-     [ 10.  26.]]
-
-Indexing base :class:`WaveletPacket2D` (as well as 1D :class:`WaveletPacket`)
-using compound path is just the same as indexing WP subnode:
-
-    >>> node = wp['a']
-    >>> print(node['a'].data)
-    [[ 10.  26.]
-     [ 10.  26.]]
-    >>> print(wp['a']['a'].data is wp['aa'].data)
-    True
-
-Following down the decomposition path:
-
-    >>> print(wp['aaa'].data)
-    [[ 36.]]
-    >>> print(wp['aaaa'].data)
-    Traceback (most recent call last):
-    ...
-    IndexError: Path length is out of range.
-
-Ups, we have reached the maximum level of decomposition for the ``'aaaa'`` path,
-which btw. was:
-
-    >>> print(wp.maxlevel)
-    3
-
-
-Now try some invalid path:
-
-    >>> print(wp['f'])
-    Traceback (most recent call last):
-    ...
-    ValueError: Subnode name must be in ['a', 'h', 'v', 'd'], not 'f'.
-
-
-Accessing Node2D's attributes:
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-:class:`WaveletPacket2D` is a tree data structure, which evaluates to a set
-of :class:`Node2D` objects. :class:`WaveletPacket2D` is just a special subclass
-of the :class:`Node2D` class (which in turn inherits from a :class:`BaseNode`,
-just like with :class:`Node` and :class:`WaveletPacket` for the 1D case.).
-
-    >>> print(wp['av'].data)
-    [[-4. -4.]
-     [-4. -4.]]
-
-    >>> print(wp['av'].path)
-    av
-
-    >>> print(wp['av'].node_name)
-    v
-
-    >>> print(wp['av'].parent.path)
-    a
-
-    >>> print(wp['av'].parent.data)
-    [[  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]]
-
-    >>> print(wp['av'].level)
-    2
-
-    >>> print(wp['av'].maxlevel)
-    3
-
-    >>> print(wp['av'].mode)
-    symmetric
-
-
-Collecting nodes
-~~~~~~~~~~~~~~~~
-
-We can get all nodes on the particular level using the
-:meth:`WaveletPacket2D.get_level` method:
-
-    * 0 level - the root `wp` node:
-
-        >>> len(wp.get_level(0))
-        1
-        >>> print([node.path for node in wp.get_level(0)])
-        ['']
-
-    * 1st level of decomposition:
-
-        >>> len(wp.get_level(1))
-        4
-        >>> print([node.path for node in wp.get_level(1)])
-        ['a', 'h', 'v', 'd']
-
-    * 2nd level of decomposition:
-
-        >>> len(wp.get_level(2))
-        16
-        >>> paths = [node.path for node in wp.get_level(2)]
-        >>> for i, path in enumerate(paths):
-        ...     if (i+1) % 4 == 0:
-        ...         print(path)
-        ...     else:
-        ...         print(path, end=' ')
-        aa ah av ad
-        ha hh hv hd
-        va vh vv vd
-        da dh dv dd
-
-    * 3rd level of decomposition:
-
-        >>> print(len(wp.get_level(3)))
-        64
-        >>> paths = [node.path for node in wp.get_level(3)]
-        >>> for i, path in enumerate(paths):
-        ...     if (i+1) % 8 == 0:
-        ...         print(path)
-        ...     else:
-        ...         print(path, end=' ')
-        aaa aah aav aad aha ahh ahv ahd
-        ava avh avv avd ada adh adv add
-        haa hah hav had hha hhh hhv hhd
-        hva hvh hvv hvd hda hdh hdv hdd
-        vaa vah vav vad vha vhh vhv vhd
-        vva vvh vvv vvd vda vdh vdv vdd
-        daa dah dav dad dha dhh dhv dhd
-        dva dvh dvv dvd dda ddh ddv ddd
-
-Note that :meth:`WaveletPacket2D.get_level` performs automatic decomposition
-until it reaches the given level.
-
-
-Reconstructing data from Wavelet Packets:
------------------------------------------
-
-Let's create a new empty 2D Wavelet Packet structure and set its nodes
-values with known data from the previous examples:
-
-    >>> new_wp = pywt.WaveletPacket2D(data=None, wavelet='db1', mode='symmetric')
-
-    >>> new_wp['vh'] = wp['vh'].data # [[0.0, 0.0], [0.0, 0.0]]
-    >>> new_wp['vv'] = wp['vh'].data # [[0.0, 0.0], [0.0, 0.0]]
-    >>> new_wp['vd'] = [[0.0, 0.0], [0.0, 0.0]]
-
-    >>> new_wp['a'] = [[3.0, 7.0, 11.0, 15.0], [3.0, 7.0, 11.0, 15.0],
-    ...                [3.0, 7.0, 11.0, 15.0], [3.0, 7.0, 11.0, 15.0]]
-    >>> new_wp['d'] = [[0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0],
-    ...                [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0]]
-
-    For convenience, :attr:`Node2D.data` gets automatically extracted from the
-    base :class:`Node2D` object:
-
-    >>> new_wp['h'] = wp['h'] # all zeros
-
-    Note: just remember to not assign to the node.data parameter directly (todo).
-
-And reconstruct the data from the ``a``, ``d``, ``vh``, ``vv``, ``vd`` and ``h``
-packets (Note that ``va`` node was not set and the WP tree is "not complete"
-- the ``va`` branch will be treated as *zero-array*):
-
-    >>> print(new_wp.reconstruct(update=False))
-    [[ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]]
-
-Now set the ``va`` node with the known values and do the reconstruction again:
-
-    >>> new_wp['va'] = wp['va'].data # [[-2.0, -2.0], [-2.0, -2.0]]
-    >>> print(new_wp.reconstruct(update=False))
-    [[ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]]
-
-which is just the same as the base sample data ``x``.
-
-Of course we can go the other way and remove nodes from the tree. If we delete
-the ``va`` node, again, we get the "not complete" tree from one of the previous
-examples:
-
-    >>> del new_wp['va']
-    >>> print(new_wp.reconstruct(update=False))
-    [[ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]
-     [ 1.5  1.5  3.5  3.5  5.5  5.5  7.5  7.5]]
-
-Just restore the node before next examples.
-
-    >>> new_wp['va'] = wp['va'].data
-
-If the ``update`` param in the :meth:`WaveletPacket2D.reconstruct` method is set
-to ``False``, the node's :attr:`Node2D.data` attribute will not be updated.
-
-    >>> print(new_wp.data)
-    None
-
-Otherwise, the :attr:`WaveletPacket2D.data` attribute will be set to the
-reconstructed value.
-
-    >>> print(new_wp.reconstruct(update=True))
-    [[ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]]
-    >>> print(new_wp.data)
-    [[ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]
-     [ 1.  2.  3.  4.  5.  6.  7.  8.]]
-
-Since we have an interesting WP structure built, it is a good occasion to
-present the :meth:`WaveletPacket2D.get_leaf_nodes` method, which collects
-non-zero leaf nodes from the WP tree:
-
-    >>> print([n.path for n in new_wp.get_leaf_nodes()])
-    ['a', 'h', 'va', 'vh', 'vv', 'vd', 'd']
-
-Passing the ``decompose = True`` parameter to the method will force the WP
-object to do a full decomposition up to the *maximum level* of decomposition:
-
-    >>> paths = [n.path for n in new_wp.get_leaf_nodes(decompose=True)]
-    >>> len(paths)
-    64
-    >>> for i, path in enumerate(paths):
-    ...     if (i+1) % 8 == 0:
-    ...         print(path)
-    ...     else:
-    ...         try:
-    ...             print(path, end=' ')
-    ...         except:
-    ...             print(path, end=' ')
-    aaa aah aav aad aha ahh ahv ahd
-    ava avh avv avd ada adh adv add
-    haa hah hav had hha hhh hhv hhd
-    hva hvh hvv hvd hda hdh hdv hdd
-    vaa vah vav vad vha vhh vhv vhd
-    vva vvh vvv vvd vda vdh vdv vdd
-    daa dah dav dad dha dhh dhv dhd
-    dva dvh dvv dvd dda ddh ddv ddd
-
-Lazy evaluation:
-----------------
-
-.. note:: This section is for demonstration of pywt internals purposes
-          only. Do not rely on the attribute access to nodes as presented in
-          this example.
-
->>> x = numpy.array([[1, 2, 3, 4, 5, 6, 7, 8]] * 8)
->>> wp = pywt.WaveletPacket2D(data=x, wavelet='db1', mode='symmetric')
-
-1) At first the wp's attribute ``a`` is ``None``
-
-   >>> print(wp.a)
-   None
-
-   **Remember that you should not rely on the attribute access.**
-
-2) During the first attempt to access the node it is computed
-   via decomposition of its parent node (the wp object itself).
-
-   >>> print(wp['a'])
-   a: [[  3.   7.  11.  15.]
-    [  3.   7.  11.  15.]
-    [  3.   7.  11.  15.]
-    [  3.   7.  11.  15.]]
-
-3) Now the ``a`` is set to the newly created node:
-
-    >>> print(wp.a)
-    a: [[  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]
-     [  3.   7.  11.  15.]]
-
-   And so is `wp.d`:
-
-    >>> print(wp.d)
-    d: [[ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]
-     [ 0.  0.  0.  0.]]
diff --git a/pyproject.toml b/pyproject.toml
index 92967a8a8..dedcbead1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -61,7 +61,7 @@ documentation = "https://pywavelets.readthedocs.io/"
 [tool.ruff]
 line-length = 88
 target-version = 'py310'
-select = [
+lint.select = [
     "I",
     "UP",
     "C4",
@@ -76,7 +76,7 @@ select = [
     "SIM118",
     "SIM2",
 ]
-ignore = [
+lint.ignore = [
      "UP038", # https://github.com/astral-sh/ruff/issues/7871
 ]
 
diff --git a/pywt/__init__.py b/pywt/__init__.py
index 3720f5504..2916bb0d1 100644
--- a/pywt/__init__.py
+++ b/pywt/__init__.py
@@ -10,8 +10,6 @@
 wavelet packets signal decomposition and reconstruction module.
 """
 
-from __future__ import division, print_function, absolute_import
-
 from ._extensions._pywt import *
 from ._functions import *
 from ._multilevel import *
diff --git a/pywt/_extensions/_pywt.pyx b/pywt/_extensions/_pywt.pyx
index 3ffe41359..ddef57500 100644
--- a/pywt/_extensions/_pywt.pyx
+++ b/pywt/_extensions/_pywt.pyx
@@ -4,7 +4,7 @@
 # See COPYING for license details.
 
 __doc__ = """Cython wrapper for low-level C wavelet transform implementation."""
-__all__ = ['MODES', 'Modes', 'DiscreteContinuousWavelet', 'Wavelet',
+__all__ = ['Modes', 'DiscreteContinuousWavelet', 'Wavelet',
            'ContinuousWavelet', 'wavelist', 'families']
 
 
diff --git a/pywt/tests/test_deprecations.py b/pywt/tests/test_deprecations.py
index 2b7d57bfb..aedaaa867 100644
--- a/pywt/tests/test_deprecations.py
+++ b/pywt/tests/test_deprecations.py
@@ -53,27 +53,6 @@ def get_mode(Modes, name):
         assert_warns(DeprecationWarning, get_mode, pywt.Modes, mode)
 
 
-def test_MODES_deprecation_new():
-    def use_MODES_new():
-        return pywt.MODES.symmetric
-
-    assert_warns(DeprecationWarning, use_MODES_new)
-
-
-def test_MODES_deprecation_old():
-    def use_MODES_old():
-        return pywt.MODES.sym
-
-    assert_warns(DeprecationWarning, use_MODES_old)
-
-
-def test_MODES_deprecation_getattr():
-    def use_MODES_new():
-        return getattr(pywt.MODES, 'symmetric')
-
-    assert_warns(DeprecationWarning, use_MODES_new)
-
-
 def test_mode_equivalence():
     old_new = [('zpd', 'zero'),
                ('cpd', 'constant'),
diff --git a/util/readthedocs/requirements.txt b/util/readthedocs/requirements.txt
index f21009667..e812cd611 100644
--- a/util/readthedocs/requirements.txt
+++ b/util/readthedocs/requirements.txt
@@ -1,13 +1,21 @@
 cython
 docutils
-jupyterlite-sphinx>=0.14.0
-jupyterlite-pyodide-kernel
+# Pyodide kernel 0.5.1 comes with Pyodide 0.27.1 and PyWavelets 1.7.0
+# see https://pyodide.org/en/stable/usage/packages-in-pyodide.html
+# and https://jupyterlite-pyodide-kernel.readthedocs.io/en/stable/#with-pyodide
+# for information on updating
+jupyterlite-pyodide-kernel==0.5.1
+jupyterlite-sphinx>=0.18.0
+jupytext
 pydata-sphinx-theme
 pytest
 matplotlib
 meson-python
+myst-nb
+nbformat
 numpy
 numpydoc
 sphinx>=7
-sphinx-togglebutton
 sphinx-copybutton
+sphinx-design
+sphinx-togglebutton