diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 22d7f47..cbea57e 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -9,12 +9,12 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macOS-latest, windows-latest]
-        python-version: [3.6, 3.7, 3.8]
+        python-version: ['3.8', '3.9', '3.10']
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
     - name: Prepare Ubuntu
diff --git a/doc/readthedocs-environment.yml b/doc/readthedocs-environment.yml
deleted file mode 100644
index d5ff057..0000000
--- a/doc/readthedocs-environment.yml
+++ /dev/null
@@ -1,14 +0,0 @@
-channels:
-  - conda-forge
-dependencies:
-  - python>=3
-  - sphinx>=1.3.6
-  - sphinx_rtd_theme
-  - sphinxcontrib-bibtex
-  - numpy
-  - scipy
-  - matplotlib>=1.5
-  - ipykernel
-  - pandoc
-  - pip:
-    - nbsphinx
diff --git a/readthedocs.yml b/readthedocs.yml
index 65bb0f4..7ff6f6a 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -1,4 +1,11 @@
-conda:
-    file: doc/readthedocs-environment.yml
+version: 2
+
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.8"
+
 python:
-    pip_install: true
+  install:
+    - requirements: requirements.txt
+    - requirements: doc/requirements.txt
diff --git a/setup.py b/setup.py
index 9a7d870..e0f0188 100644
--- a/setup.py
+++ b/setup.py
@@ -13,7 +13,7 @@
     version=__version__,
     packages=find_packages(),
     install_requires=[
-        'numpy!=1.11.0',  # https://github.com/sfstoolbox/sfs-python/issues/11
+        'numpy',
         'scipy',
     ],
     author="SFS Toolbox Developers",
@@ -31,8 +31,9 @@
         "Operating System :: OS Independent",
         "Programming Language :: Python",
         "Programming Language :: Python :: 3",
-        "Programming Language :: Python :: 3.6",
-        "Programming Language :: Python :: 3.7",
+        "Programming Language :: Python :: 3.8",
+        "Programming Language :: Python :: 3.9",
+        "Programming Language :: Python :: 3.10",
         "Programming Language :: Python :: 3 :: Only",
         "Topic :: Scientific/Engineering",
     ],
diff --git a/sfs/fd/wfs.py b/sfs/fd/wfs.py
index 44fb2c6..840b59e 100644
--- a/sfs/fd/wfs.py
+++ b/sfs/fd/wfs.py
@@ -32,7 +32,6 @@ def plot(d, selection, secondary_source):
 
 """
 import numpy as _np
-from numpy.core.umath_tests import inner1d as _inner1d
 from scipy.special import hankel2 as _hankel2
 
 from . import secondary_source_line as _secondary_source_line
@@ -91,7 +90,7 @@ def line_2d(omega, x0, n0, xs, *, c=None):
     k = _util.wavenumber(omega, c)
     ds = x0 - xs
     r = _np.linalg.norm(ds, axis=1)
-    d = -1j/2 * k * _inner1d(ds, n0) / r * _hankel2(1, k * r)
+    d = -1j/2 * k * _util._inner1d(ds, n0) / r * _hankel2(1, k * r)
     selection = _util.source_selection_line(n0, x0, xs)
     return d, selection, _secondary_source_line(omega, c)
 
@@ -147,7 +146,7 @@ def _point(omega, x0, n0, xs, *, c=None):
     k = _util.wavenumber(omega, c)
     ds = x0 - xs
     r = _np.linalg.norm(ds, axis=1)
-    d = 1j * k * _inner1d(ds, n0) / r ** (3 / 2) * _np.exp(-1j * k * r)
+    d = 1j * k * _util._inner1d(ds, n0) / r ** (3 / 2) * _np.exp(-1j * k * r)
     selection = _util.source_selection_point(n0, x0, xs)
     return d, selection, _secondary_source_point(omega, c)
 
@@ -234,7 +233,7 @@ def point_25d(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None):
         preeq_25d(omega, omalias, c) *
         _np.sqrt(8 * _np.pi) *
         _np.sqrt((r * s) / (r + s)) *
-        _inner1d(n0, ds) / s *
+        _util._inner1d(n0, ds) / s *
         _np.exp(-1j * k * s) / (4 * _np.pi * s))
     selection = _util.source_selection_point(n0, x0, xs)
     return d, selection, _secondary_source_point(omega, c)
@@ -316,7 +315,7 @@ def point_25d_legacy(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None):
     r = _np.linalg.norm(ds, axis=1)
     d = (
         preeq_25d(omega, omalias, c) *
-        _np.sqrt(_np.linalg.norm(xref - x0)) * _inner1d(ds, n0) /
+        _np.sqrt(_np.linalg.norm(xref - x0)) * _util._inner1d(ds, n0) /
         r ** (3 / 2) * _np.exp(-1j * k * r))
     selection = _util.source_selection_point(n0, x0, xs)
     return d, selection, _secondary_source_point(omega, c)
@@ -499,7 +498,7 @@ def _focused(omega, x0, n0, xs, ns, *, c=None):
     k = _util.wavenumber(omega, c)
     ds = x0 - xs
     r = _np.linalg.norm(ds, axis=1)
-    d = 1j * k * _inner1d(ds, n0) / r ** (3 / 2) * _np.exp(1j * k * r)
+    d = 1j * k * _util._inner1d(ds, n0) / r ** (3 / 2) * _np.exp(1j * k * r)
     selection = _util.source_selection_focused(ns, x0, xs)
     return d, selection, _secondary_source_point(omega, c)
 
@@ -569,7 +568,7 @@ def focused_25d(omega, x0, n0, xs, ns, *, xref=[0, 0, 0], c=None,
     r = _np.linalg.norm(ds, axis=1)
     d = (
         preeq_25d(omega, omalias, c) *
-        _np.sqrt(_np.linalg.norm(xref - x0)) * _inner1d(ds, n0) /
+        _np.sqrt(_np.linalg.norm(xref - x0)) * _util._inner1d(ds, n0) /
         r ** (3 / 2) * _np.exp(1j * k * r))
     selection = _util.source_selection_focused(ns, x0, xs)
     return d, selection, _secondary_source_point(omega, c)
diff --git a/sfs/plot2d.py b/sfs/plot2d.py
index efadb9c..63f6996 100644
--- a/sfs/plot2d.py
+++ b/sfs/plot2d.py
@@ -18,13 +18,17 @@ def _register_cmap_clip(name, original_cmap, alpha):
         cmap = LinearSegmentedColormap.from_list(name, cdata)
     cmap.set_over([alpha * c + 1 - alpha for c in cmap(1.0)[:3]])
     cmap.set_under([alpha * c + 1 - alpha for c in cmap(0.0)[:3]])
-    _plt.cm.register_cmap(cmap=cmap)
+    _plt.colormaps.register(cmap=cmap)
 
 
 # The 'coolwarm' colormap is based on the paper
 # "Diverging Color Maps for Scientific Visualization" by Kenneth Moreland
 # http://www.sandia.gov/~kmorel/documents/ColorMaps/
-_register_cmap_clip('coolwarm_clip', 'coolwarm', 0.7)
+# already registered in MPL 3.9.0
+try:
+    _register_cmap_clip('coolwarm_clip', 'coolwarm', 0.7)
+except ImportError:
+    pass
 
 
 def _register_cmap_transparent(name, color):
@@ -36,7 +40,7 @@ def _register_cmap_transparent(name, color):
              'blue': ((0, blue, blue), (1, blue, blue)),
              'alpha': ((0, 0, 0), (1, 1, 1))}
     cmap = LinearSegmentedColormap(name, cdict)
-    _plt.cm.register_cmap(cmap=cmap)
+    _plt.colormaps.register(cmap=cmap)
 
 
 _register_cmap_transparent('blacktransparent', 'black')
@@ -285,8 +289,8 @@ def amplitude(p, grid, *, xnorm=None, cmap='coolwarm_clip',
     elif plotting_plane == 'yz':
         x, y = grid[[1, 2]]
 
-    dx = 0.5 * x.ptp() / p.shape[0]
-    dy = 0.5 * y.ptp() / p.shape[1]
+    dx = 0.5 * _np.ptp(x) / p.shape[0]
+    dy = 0.5 * _np.ptp(y) / p.shape[1]
 
     if ax is None:
         ax = _plt.gca()
diff --git a/sfs/td/wfs.py b/sfs/td/wfs.py
index 3b59301..3e11208 100644
--- a/sfs/td/wfs.py
+++ b/sfs/td/wfs.py
@@ -44,7 +44,6 @@ def plot(d, selection, secondary_source, t=0):
 
 """
 import numpy as _np
-from numpy.core.umath_tests import inner1d as _inner1d
 
 from . import apply_delays as _apply_delays
 from . import secondary_source_point as _secondary_source_point
@@ -119,8 +118,8 @@ def plane_25d(x0, n0, n=[0, 1, 0], xref=[0, 0, 0], c=None):
     n = _util.normalize_vector(n)
     xref = _util.asarray_1d(xref)
     g0 = _np.sqrt(2 * _np.pi * _np.linalg.norm(xref - x0, axis=1))
-    delays = _inner1d(n, x0) / c
-    weights = 2 * g0 * _inner1d(n, n0)
+    delays = _util._inner1d(n, x0) / c
+    weights = 2 * g0 * _util._inner1d(n, n0)
     selection = _util.source_selection_plane(n0, n)
     return delays, weights, selection, _secondary_source_point(c)
 
@@ -208,7 +207,7 @@ def point_25d(x0, n0, xs, xref=[0, 0, 0], c=None):
     g0 *= _np.sqrt((x0xs_n*x0xref_n)/(x0xs_n+x0xref_n))
 
     delays = x0xs_n/c
-    weights = g0*_inner1d(x0xs, n0)
+    weights = g0 * _util._inner1d(x0xs, n0)
     selection = _util.source_selection_point(n0, x0, xs)
     return delays, weights, selection, _secondary_source_point(c)
 
@@ -296,7 +295,7 @@ def point_25d_legacy(x0, n0, xs, xref=[0, 0, 0], c=None):
     ds = x0 - xs
     r = _np.linalg.norm(ds, axis=1)
     delays = r/c
-    weights = g0 * _inner1d(ds, n0) / (2 * _np.pi * r**(3/2))
+    weights = g0 * _util._inner1d(ds, n0) / (2 * _np.pi * r**(3/2))
     selection = _util.source_selection_point(n0, x0, xs)
     return delays, weights, selection, _secondary_source_point(c)
 
@@ -379,7 +378,7 @@ def focused_25d(x0, n0, xs, ns, xref=[0, 0, 0], c=None):
     g0 = _np.sqrt(_np.linalg.norm(xref - x0, axis=1)
                   / (_np.linalg.norm(xref - x0, axis=1) + r))
     delays = -r/c
-    weights = g0 * _inner1d(ds, n0) / (2 * _np.pi * r**(3/2))
+    weights = g0 * _util._inner1d(ds, n0) / (2 * _np.pi * r**(3/2))
     selection = _util.source_selection_focused(ns, x0, xs)
     return delays, weights, selection, _secondary_source_point(c)
 
diff --git a/sfs/util.py b/sfs/util.py
index c15358f..de2ead5 100644
--- a/sfs/util.py
+++ b/sfs/util.py
@@ -6,7 +6,6 @@
 
 import collections
 import numpy as np
-from numpy.core.umath_tests import inner1d
 from scipy.special import spherical_jn, spherical_yn
 from . import default
 
@@ -576,7 +575,7 @@ def source_selection_point(n0, x0, xs):
     x0 = asarray_of_rows(x0)
     xs = asarray_1d(xs)
     ds = x0 - xs
-    return inner1d(ds, n0) >= default.selection_tolerance
+    return _inner1d(ds, n0) >= default.selection_tolerance
 
 
 def source_selection_line(n0, x0, xs):
@@ -598,7 +597,7 @@ def source_selection_focused(ns, x0, xs):
     xs = asarray_1d(xs)
     ns = normalize_vector(ns)
     ds = xs - x0
-    return inner1d(ns, ds) >= default.selection_tolerance
+    return _inner1d(ns, ds) >= default.selection_tolerance
 
 
 def source_selection_all(N):
@@ -646,3 +645,8 @@ def max_order_spherical_harmonics(N):
 
     """
     return int(np.sqrt(N) - 1)
+
+
+def _inner1d(arr1, arr2):
+    # https://github.com/numpy/numpy/issues/10815#issuecomment-376847774
+    return (arr1 * arr2).sum(axis=1)