diff --git a/.gitignore b/.gitignore
index cebac22c5..f862d8a3a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,6 +16,7 @@ venv/
 ENV/
 env.bak/
 venv.bak/
+env8
 
 # automatically generated by setuptools-scm
 pysindy/version.py
diff --git a/pyproject.toml b/pyproject.toml
index 1e7f22cf1..9d614262e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -50,7 +50,7 @@ docs = [
     "ipython",
     "pandoc",
     "sphinx-rtd-theme",
-    "sphinx==5.3.0",
+    "sphinx==7.1.2",
     "sphinxcontrib-apidoc",
     "nbsphinx"
 ]
diff --git a/pysindy/differentiation/finite_difference.py b/pysindy/differentiation/finite_difference.py
index 7516b44e9..0c44917e9 100644
--- a/pysindy/differentiation/finite_difference.py
+++ b/pysindy/differentiation/finite_difference.py
@@ -1,7 +1,12 @@
+from typing import List
+from typing import Union
+
 import numpy as np
+from numpy.typing import NDArray
 from scipy.special import factorial
 
 from .base import BaseDifferentiation
+from pysindy.utils.axes import AxesArray
 
 
 class FiniteDifference(BaseDifferentiation):
@@ -55,7 +60,7 @@ class FiniteDifference(BaseDifferentiation):
     def __init__(
         self,
         order=2,
-        d=1,
+        d: int = 1,
         axis=0,
         is_uniform=False,
         drop_endpoints=False,
@@ -85,22 +90,27 @@ def __init__(
 
     def _coefficients(self, t):
         nt = len(t)
-        self.stencil_inds = np.array(
-            [np.arange(i, nt - self.n_stencil + i + 1) for i in range(self.n_stencil)]
+        self.stencil_inds = AxesArray(
+            np.array(
+                [
+                    np.arange(i, nt - self.n_stencil + i + 1)
+                    for i in range(self.n_stencil)
+                ]
+            ),
+            {"ax_offset": 0, "ax_ti": 1},
+        )
+        self.stencil = AxesArray(
+            np.transpose(t[self.stencil_inds]), {"ax_time": 0, "ax_offset": 1}
         )
-        self.stencil = np.transpose(t[self.stencil_inds])
-
         pows = np.arange(self.n_stencil)[np.newaxis, :, np.newaxis]
-        matrices = (
+        dt_endpoints = (
             self.stencil
-            - t[
-                (self.n_stencil - 1) // 2 : -(self.n_stencil - 1) // 2,
-                np.newaxis,
-            ]
-        )[:, np.newaxis, :] ** pows
-        b = np.zeros(self.n_stencil)
-        b[self.d] = factorial(self.d)
-        return np.linalg.solve(matrices, [b])
+            - t[(self.n_stencil - 1) // 2 : -(self.n_stencil - 1) // 2, "offset"]
+        )
+        matrices = dt_endpoints[:, "power", :] ** pows
+        b = AxesArray(np.zeros((1, self.n_stencil)), {"ax_time": 0, "ax_power": 1})
+        b[0, self.d] = factorial(self.d)
+        return np.linalg.solve(matrices, b)
 
     def _coefficients_boundary_forward(self, t):
         # use the same stencil for each boundary point,
@@ -202,23 +212,30 @@ def _constant_coefficients(self, dt):
 
     def _accumulate(self, coeffs, x):
         # slice to select the stencil indices
-        s = [slice(None)] * len(x.shape)
+        s = [slice(None)] * x.ndim
         s[self.axis] = self.stencil_inds
 
-        # a new axis is introduced after self.axis for the stencil indices
+        # a new axis is introduced before self.axis for the stencil indices
         # To contract with the coefficients, roll by -self.axis to put it first
         # Then roll back by self.axis to return the order
-        trans = np.roll(np.arange(len(x.shape) + 1), -self.axis)
+        trans = np.roll(np.arange(x.ndim + 1), -self.axis)
+        # TODO: assign x's axes much earlier in the call stack
+        x = AxesArray(x, {"ax_unk": list(range(x.ndim))})
+        x_expanded = AxesArray(
+            np.transpose(x[tuple(s)], axes=trans), x.insert_axis(0, "ax_offset")
+        )
         return np.transpose(
             np.einsum(
                 "ij...,ij->j...",
-                np.transpose(x[tuple(s)], axes=trans),
+                x_expanded,
                 np.transpose(coeffs),
             ),
-            np.roll(np.arange(len(x.shape)), self.axis),
+            np.roll(np.arange(x.ndim), self.axis),
         )
 
-    def _differentiate(self, x, t):
+    def _differentiate(
+        self, x: NDArray, t: Union[NDArray, float, List[float]]
+    ) -> NDArray:
         """
         Apply finite difference method.
         """
@@ -249,6 +266,7 @@ def _differentiate(self, x, t):
                     s[self.axis] = slice(start, stop)
                     interior = interior + x[tuple(s)] * coeffs[i]
         else:
+            t = AxesArray(np.array(t), axes={"ax_time": 0})
             coeffs = self._coefficients(t)
             interior = self._accumulate(coeffs, x)
         s[self.axis] = slice((self.n_stencil - 1) // 2, -(self.n_stencil - 1) // 2)
diff --git a/pysindy/feature_library/base.py b/pysindy/feature_library/base.py
index 16149b27c..54697da45 100644
--- a/pysindy/feature_library/base.py
+++ b/pysindy/feature_library/base.py
@@ -63,10 +63,9 @@ def correct_shape(self, x: AxesArray):
         return x
 
     def calc_trajectory(self, diff_method, x, t):
-        axes = x.__dict__
         x_dot = diff_method(x, t=t)
-        x = AxesArray(diff_method.smoothed_x_, axes)
-        return x, AxesArray(x_dot, axes)
+        x = AxesArray(diff_method.smoothed_x_, x.axes)
+        return x, AxesArray(x_dot, x.axes)
 
     def get_spatial_grid(self):
         return None
@@ -337,7 +336,7 @@ def __init__(
         self.libraries = libraries
         self.inputs_per_library = inputs_per_library
 
-    def _combinations(self, lib_i, lib_j):
+    def _combinations(self, lib_i: AxesArray, lib_j: AxesArray) -> AxesArray:
         """
         Compute combinations of the numerical libraries.
 
@@ -351,7 +350,7 @@ def _combinations(self, lib_i, lib_j):
             lib_i.shape[lib_i.ax_coord] * lib_j.shape[lib_j.ax_coord]
         )
         lib_full = np.reshape(
-            lib_i[..., :, np.newaxis] * lib_j[..., np.newaxis, :],
+            lib_i[..., :, "coord"] * lib_j[..., "coord", :],
             shape,
         )
 
diff --git a/pysindy/feature_library/generalized_library.py b/pysindy/feature_library/generalized_library.py
index 3e5e24055..29834c2a8 100644
--- a/pysindy/feature_library/generalized_library.py
+++ b/pysindy/feature_library/generalized_library.py
@@ -237,7 +237,7 @@ def transform(self, x_full):
                 else:
                     xps.append(lib.transform([x])[0])
 
-            xp = AxesArray(np.concatenate(xps, axis=xps[0].ax_coord), xps[0].__dict__)
+            xp = AxesArray(np.concatenate(xps, axis=xps[0].ax_coord), xps[0].axes)
             xp_full = xp_full + [xp]
         return xp_full
 
diff --git a/pysindy/feature_library/pde_library.py b/pysindy/feature_library/pde_library.py
index cc9ec0b8c..bce1b1a48 100644
--- a/pysindy/feature_library/pde_library.py
+++ b/pysindy/feature_library/pde_library.py
@@ -234,13 +234,7 @@ def get_feature_names(self, input_features=None):
         def derivative_string(multiindex):
             ret = ""
             for axis in range(self.ind_range):
-                if self.implicit_terms and (
-                    axis
-                    in [
-                        self.spatiotemporal_grid.ax_time,
-                        self.spatiotemporal_grid.ax_sample,
-                    ]
-                ):
+                if self.implicit_terms and (axis == self.spatiotemporal_grid.ax_time,):
                     str_deriv = "t"
                 else:
                     str_deriv = str(axis + 1)
@@ -345,7 +339,7 @@ def transform(self, x_full):
 
             # derivative terms
             shape[-1] = n_features * self.num_derivatives
-            library_derivatives = np.empty(shape, dtype=x.dtype)
+            library_derivatives = AxesArray(np.empty(shape, dtype=x.dtype), x.axes)
             library_idx = 0
             for multiindex in self.multiindices:
                 derivs = x
@@ -395,8 +389,8 @@ def transform(self, x_full):
                     library_idx : library_idx
                     + n_library_terms * self.num_derivatives * n_features,
                 ] = np.reshape(
-                    library_functions[..., np.newaxis, :]
-                    * library_derivatives[..., :, np.newaxis],
+                    library_functions[..., "coord", :]
+                    * library_derivatives[..., :, "coord"],
                     shape,
                 )
                 library_idx += n_library_terms * self.num_derivatives * n_features
diff --git a/pysindy/feature_library/polynomial_library.py b/pysindy/feature_library/polynomial_library.py
index 75dbf5637..e62af38bd 100644
--- a/pysindy/feature_library/polynomial_library.py
+++ b/pysindy/feature_library/polynomial_library.py
@@ -225,7 +225,7 @@ def transform(self, x_full):
                         dtype=x.dtype,
                         order=self.order,
                     ),
-                    x.__dict__,
+                    x.axes,
                 )
                 for i, comb in enumerate(combinations):
                     xp[..., i] = x[..., comb].prod(-1)
diff --git a/pysindy/feature_library/sindy_pi_library.py b/pysindy/feature_library/sindy_pi_library.py
index 8d5f054a7..f45cf567f 100644
--- a/pysindy/feature_library/sindy_pi_library.py
+++ b/pysindy/feature_library/sindy_pi_library.py
@@ -404,5 +404,5 @@ def transform(self, x_full):
                                     *[x[:, comb] for comb in f_combs]
                                 ) * f_dot(*[x_dot[:, comb] for comb in f_dot_combs])
                                 library_idx += 1
-            xp_full = xp_full + [AxesArray(xp, x.__dict__)]
+            xp_full = xp_full + [AxesArray(xp, x.axes)]
         return xp_full
diff --git a/pysindy/feature_library/weak_pde_library.py b/pysindy/feature_library/weak_pde_library.py
index 65b551da4..0f41ede84 100644
--- a/pysindy/feature_library/weak_pde_library.py
+++ b/pysindy/feature_library/weak_pde_library.py
@@ -8,6 +8,7 @@
 from sklearn.utils.validation import check_is_fitted
 
 from ..utils import AxesArray
+from ..utils import comprehend_axes
 from .base import BaseFeatureLibrary
 from .base import x_sequence_or_item
 from .polynomial_library import PolynomialLibrary
@@ -218,7 +219,10 @@ def __init__(
 
         self.num_derivatives = num_derivatives
         self.multiindices = multiindices
-        self.spatiotemporal_grid = spatiotemporal_grid
+
+        self.spatiotemporal_grid = AxesArray(
+            spatiotemporal_grid, axes=comprehend_axes(spatiotemporal_grid)
+        )
 
         # Weak form checks and setup
         self._weak_form_setup()
@@ -228,12 +232,14 @@ def _weak_form_setup(self):
         L_xt = xt2 - xt1
         if self.H_xt is not None:
             if np.isscalar(self.H_xt):
-                self.H_xt = np.array(self.grid_ndim * [self.H_xt])
+                self.H_xt = AxesArray(
+                    np.array(self.grid_ndim * [self.H_xt]), {"ax_coord": 0}
+                )
             if self.grid_ndim != len(self.H_xt):
                 raise ValueError(
                     "The user-defined grid (spatiotemporal_grid) and "
                     "the user-defined sizes of the subdomains for the "
-                    "weak form, do not have the same # of spatiotemporal "
+                    "weak form do not have the same # of spatiotemporal "
                     "dimensions. For instance, if spatiotemporal_grid is 4D, "
                     "then H_xt should be a 4D list of the subdomain lengths."
                 )
@@ -258,8 +264,8 @@ def _weak_form_setup(self):
         self._set_up_weights()
 
     def _get_spatial_endpoints(self):
-        x1 = np.zeros(self.grid_ndim)
-        x2 = np.zeros(self.grid_ndim)
+        x1 = AxesArray(np.zeros(self.grid_ndim), {"ax_coord": 0})
+        x2 = AxesArray(np.zeros(self.grid_ndim), {"ax_coord": 0})
         for i in range(self.grid_ndim):
             inds = [slice(None)] * (self.grid_ndim + 1)
             for j in range(self.grid_ndim):
@@ -279,7 +285,9 @@ def _set_up_weights(self):
 
         # Sample the random domain centers
         xt1, xt2 = self._get_spatial_endpoints()
-        domain_centers = np.zeros((self.K, self.grid_ndim))
+        domain_centers = AxesArray(
+            np.zeros((self.K, self.grid_ndim)), {"ax_sample": 0, "ax_coord": 1}
+        )
         for i in range(self.grid_ndim):
             domain_centers[:, i] = np.random.uniform(
                 xt1[i] + self.H_xt[i], xt2[i] - self.H_xt[i], size=self.K
@@ -294,15 +302,12 @@ def _set_up_weights(self):
                 s = [0] * (self.grid_ndim + 1)
                 s[i] = slice(None)
                 s[-1] = i
-                newinds = np.intersect1d(
-                    np.where(
-                        self.spatiotemporal_grid[tuple(s)]
-                        >= domain_centers[k][i] - self.H_xt[i]
-                    ),
-                    np.where(
-                        self.spatiotemporal_grid[tuple(s)]
-                        <= domain_centers[k][i] + self.H_xt[i]
-                    ),
+                ax_vals = self.spatiotemporal_grid[tuple(s)]
+                cell_left = domain_centers[k][i] - self.H_xt[i]
+                cell_right = domain_centers[k][i] + self.H_xt[i]
+                newinds = AxesArray(
+                    ((ax_vals > cell_left) & (ax_vals < cell_right)).nonzero()[0],
+                    ax_vals.axes,
                 )
                 # If less than two indices along any axis, resample
                 if len(newinds) < 2:
@@ -319,6 +324,7 @@ def _set_up_weights(self):
                 self.inds_k = self.inds_k + [inds]
                 k = k + 1
 
+        # TODO: fix meaning of axes in XT_k
         # Values of the spatiotemporal grid on the domain cells
         XT_k = [
             self.spatiotemporal_grid[np.ix_(*self.inds_k[k])] for k in range(self.K)
@@ -441,6 +447,11 @@ def _set_up_weights(self):
                 )
             weights1 = weights1 + [weights2]
 
+        # TODO: get rest of code to work with AxesArray.  Too unsure of
+        # which axis labels to use at this point to continue
+        tweights = [np.asarray(arr) for arr in tweights]
+        weights0 = [np.asarray(arr) for arr in weights0]
+        weights1 = [[np.asarray(arr) for arr in sublist] for sublist in weights1]
         # Product weights over the axes for time derivatives, shaped as inds_k
         self.fulltweights = []
         deriv = np.zeros(self.grid_ndim)
diff --git a/pysindy/optimizers/base.py b/pysindy/optimizers/base.py
index 45d4842b2..614341b54 100644
--- a/pysindy/optimizers/base.py
+++ b/pysindy/optimizers/base.py
@@ -144,7 +144,8 @@ def fit(self, x_, y, sample_weight=None, **reduce_kws):
         self : returns an instance of self
         """
         x_ = AxesArray(np.asarray(x_), {"ax_sample": 0, "ax_coord": 1})
-        y = AxesArray(np.asarray(y), {"ax_sample": 0, "ax_coord": 1})
+        y_axes = {"ax_sample": 0} if y.ndim == 1 else {"ax_sample": 0, "ax_coord": 1}
+        y = AxesArray(np.asarray(y), y_axes)
         x_, y = drop_nan_samples(x_, y)
         x_, y = check_X_y(x_, y, accept_sparse=[], y_numeric=True, multi_output=True)
 
diff --git a/pysindy/utils/axes.py b/pysindy/utils/axes.py
index ad0f79040..ed27957aa 100644
--- a/pysindy/utils/axes.py
+++ b/pysindy/utils/axes.py
@@ -1,63 +1,357 @@
+"""
+A module that defines one external class, AxesArray, to act like a numpy array
+but keep track of axis definitions.  It aims to allow meaningful replacement
+of magic numbers for axis conventions in code.  E.g::
+
+   import numpy as np
+
+   arr = AxesArray(np.ones((2,3,4)), {"ax_time": 0, "ax_spatial": [1, 2]})
+   print(arr.axes)
+   print(arr.ax_time)
+   print(arr.n_time)
+   print(arr.ax_spatial)
+   print(arr.n_spatial)
+
+Would show::
+
+   {"ax_time": 0, "ax_spatial": [1, 2]}
+   0
+   2
+   [1, 2]
+   [3, 4]
+
+It is up to the user to handle the ``list[int] | int`` return values, but this
+module has several functions to deal with the axes dictionary, internally
+referred to as type ``CompatDict[T]``:
+
+Appending an item to a ``CompatDict[T]``
+   :py:func:`compat_dict_append`
+
+Generating a ``CompatDict[int]`` of axes from list of axes names:
+   :py:func:`fwd_from_names`
+
+Create new ``CompatDict[int]`` from this ``AxesArray`` with new axis/axes added:
+   :py:meth:`AxesArray.insert_axis`
+
+Create new ``CompatDict[int]`` from this ``AxesArray`` with axis/axes removed:
+   :py:meth:`AxesArray.remove_axis`
+
+
+.. todo::
+
+   Add developer documentation here.
+
+The recommended way to refactor existing code to use AxesArrays is to add them
+at the lowest level possible.  Enter debug mode and see how long the expected
+axes persist throughout array operations.  When AxesArray loses track of the
+correct axes, re-assign them with an AxesArray constructor (which only uses a
+view of the data).
+
+Starting at the macro level runs the risk of triggering a great deal of errors
+from unimplemented functions.
+"""
+from __future__ import annotations
+
+import copy
+import warnings
+from enum import Enum
+from typing import Collection
+from typing import Dict
+from typing import get_args
 from typing import List
+from typing import Literal
+from typing import NewType
+from typing import Optional
+from typing import Sequence
+from typing import Tuple
+from typing import TypeVar
+from typing import Union
 
 import numpy as np
+from numpy.typing import NDArray
 from sklearn.base import TransformerMixin
 
 HANDLED_FUNCTIONS = {}
 
+AxesWarning = type("AxesWarning", (SyntaxWarning,), {})
+BasicIndexer = Union[slice, int, type(Ellipsis), None, str]
+Indexer = Union[BasicIndexer, NDArray, List]
+StandardIndexer = Union[slice, int, None, NDArray[np.dtype(int)]]
+OldIndex = NewType("OldIndex", int)  # Before moving advanced axes adajent
+KeyIndex = NewType("KeyIndex", int)
+NewIndex = NewType("NewIndex", int)
+T = TypeVar("T", bound=int)  # TODO: Bind to a non-sequence after type-negation PEP
+ItemOrList = Union[T, List[T]]
+CompatDict = Dict[str, ItemOrList[T]]
+
+
+class _Sentinels(Enum):
+    ADV_NAME = object()
+    ADV_REMOVE = object()
+
+
+class _AxisMapping:
+    """Convenience wrapper for a two-way map between axis names and indexes."""
+
+    fwd_map: Dict[str, List[int]]
+    reverse_map: Dict[int, str]
+
+    def __init__(
+        self,
+        axes: dict[str, Union[int, Sequence[int]]],
+        in_ndim: int,
+    ):
+        self.fwd_map = {}
+        self.reverse_map = {}
+
+        def coerce_sequence(obj):
+            if isinstance(obj, Sequence):
+                return sorted(obj)
+            return [obj]
+
+        for ax_name, ax_ids in axes.items():
+            ax_ids = coerce_sequence(ax_ids)
+            self.fwd_map[ax_name] = ax_ids
+            for ax_id in ax_ids:
+                old_name = self.reverse_map.get(ax_id)
+                if old_name is not None:
+                    raise ValueError(f"Assigned multiple definitions to axis {ax_id}")
+                if ax_id >= in_ndim:
+                    raise ValueError(
+                        f"Assigned definition to axis {ax_id}, but array only has"
+                        f" {in_ndim} axes"
+                    )
+                self.reverse_map[ax_id] = ax_name
+        if len(self.reverse_map) != in_ndim:
+            warnings.warn(
+                f"{len(self.reverse_map)} axes labeled for array with {in_ndim} axes",
+                AxesWarning,
+            )
+
+    @staticmethod
+    def _compat_axes(in_dict: Dict[str, List[int]]) -> Dict[str, Union[list[int], int]]:
+        """Like fwd_map, but unpack single-element axis lists"""
+        axes = {}
+        for k, v in in_dict.items():
+            if len(v) == 1:
+                axes[k] = v[0]
+            else:
+                axes[k] = v
+        return axes
+
+    @property
+    def compat_axes(self):
+        return self._compat_axes(self.fwd_map)
+
+    def remove_axis(self, axis: Union[Collection[int], int, None] = None):
+        """Create an axes dict from self with specified axis or axes
+        removed and all greater axes decremented.  This can be passed to
+        the constructor to create a new _AxisMapping
+
+        Arguments:
+            axis: the axis index or axes indexes to remove.  By numpy
+            ufunc convention, axis=None (default) removes _all_ axes.
+        """
+        if axis is None:
+            return {}
+        new_axes = copy.deepcopy(self.fwd_map)
+        in_ndim = len(self.reverse_map)
+        if not isinstance(axis, Collection):
+            axis = [axis]
+        axis = [ax_id if ax_id >= 0 else (self.ndim + ax_id) for ax_id in axis]
+        for cum_shift, orig_ax_remove in enumerate(sorted(axis)):
+            remove_ax_name = self.reverse_map[orig_ax_remove]
+            curr_ax_remove = orig_ax_remove - cum_shift
+            if len(new_axes[remove_ax_name]) == 1:
+                new_axes.pop(remove_ax_name)
+            else:
+                new_axes[remove_ax_name].remove(curr_ax_remove)
+            for old_ax_dec in range(curr_ax_remove + 1, in_ndim - cum_shift):
+                orig_ax_dec = old_ax_dec + cum_shift
+                ax_dec_name = self.reverse_map[orig_ax_dec]
+                new_axes[ax_dec_name].remove(old_ax_dec)
+                new_axes[ax_dec_name].append(old_ax_dec - 1)
+        return self._compat_axes(new_axes)
+
+    def insert_axis(self, axis: Union[Collection[int], int], new_name: str):
+        """Create an axes dict from self with specified axis or axes
+        added and all greater axes incremented.
+
+        Arguments:
+            axis: the axis index or axes indexes to add.
+
+        Todo:
+            May be more efficient to determine final axis-to-axis
+            mapping, then apply, rather than apply changes after each
+            axis insert.
+        """
+        new_axes = copy.deepcopy(self.fwd_map)
+        in_ndim = len(self.reverse_map)
+        if not isinstance(axis, Collection):
+            axis = [axis]
+        for cum_shift, ax in enumerate(sorted(axis)):
+            if new_name in new_axes.keys():
+                new_axes[new_name].append(ax)
+            else:
+                new_axes[new_name] = [ax]
+            for ax_id in range(ax, in_ndim + cum_shift):
+                ax_name = self.reverse_map[ax_id - cum_shift]
+                new_axes[ax_name].remove(ax_id)
+                new_axes[ax_name].append(ax_id + 1)
+        return self._compat_axes(new_axes)
+
+    @property
+    def ndim(self):
+        return len(self.reverse_map)
+
 
 class AxesArray(np.lib.mixins.NDArrayOperatorsMixin, np.ndarray):
     """A numpy-like array that keeps track of the meaning of its axes.
 
+    Limitations:
+
+    * Not all numpy functions, such as ``np.flatten()``, have an
+      implementation for ``AxesArray``.  In such cases a regular numpy array
+      is returned.
+    * For functions that are implemented for `AxesArray`, such as
+      ``np.reshape()``, use the numpy function rather than the bound
+      method (e.g. ``arr.reshape``)
+    * Such functions may raise ``ValueError`` where numpy would not, when
+      it is impossible to determine the output axis labels.
+
+    Current array function implementations:
+
+    * ``np.concatenate``
+    * ``np.reshape``
+    * ``np.transpose``
+    * ``np.linalg.solve``
+    * ``np.einsum``
+    * ``np.tensordot``
+
+    Indexing:
+        AxesArray supports all of the basic and advanced indexing of numpy
+        arrays, with the addition that new axes can be inserted with a string
+        name for the axis.  E.g. ``arr = arr[..., "lineno"]`` will add a
+        length-one axis at the end, along with the properties ``arr.ax_lineno``
+        and ``arr.n_lineno``.  If ``None`` or ``np.newaxis`` are passed, the
+        axis is named "unk".
+
     Parameters:
-        input_array (array-like): the data to create the array.
-        axes (dict): A dictionary of axis labels to shape indices.
-            Allowed keys:
-            -  ax_time: int
-            -  ax_coord: int
-            -  ax_sample: int
-            -  ax_spatial: List[int]
+        input_array: the data to create the array.
+        axes: A dictionary of axis labels to shape indices.  Axes labels must
+            be of the format "ax_name".  indices can be either an int or a
+            list of ints.
+
+    Attributes:
+        axes: dictionary of axis name to dimension index/indices
+        ax_<ax_name>: lookup ax_name in axes
+        n_<ax_name>: lookup shape of subarray defined by ax_name
 
     Raises:
-        AxesWarning if axes does not match shape of input_array
+        AxesWarning if axes does not match shape of input_array.
+        ValueError if assigning the same axis index to multiple meanings or
+            assigning an axis beyond ndim.
+
     """
 
-    def __new__(cls, input_array, axes):
+    _ax_map: _AxisMapping
+
+    def __new__(cls, input_array: NDArray, axes: CompatDict[int]):
         obj = np.asarray(input_array).view(cls)
-        defaults = {
-            "ax_time": None,
-            "ax_coord": None,
-            "ax_sample": None,
-            "ax_spatial": [],
-        }
-        if axes is None:
-            return obj
-        obj.__dict__.update({**defaults, **axes})
+        in_ndim = len(input_array.shape)
+        obj._ax_map = _AxisMapping(axes, in_ndim)
         return obj
 
-    def __array_finalize__(self, obj) -> None:
-        if obj is None:
-            return
-        self.ax_time = getattr(obj, "ax_time", None)
-        self.ax_coord = getattr(obj, "ax_coord", None)
-        self.ax_sample = getattr(obj, "ax_sample", None)
-        self.ax_spatial = getattr(obj, "ax_spatial", [])
-
     @property
-    def n_spatial(self):
-        return tuple(self.shape[ax] for ax in self.ax_spatial)
+    def axes(self):
+        return self._ax_map.compat_axes
 
     @property
-    def n_time(self):
-        return self.shape[self.ax_time] if self.ax_time is not None else 1
+    def _reverse_map(self):
+        return self._ax_map.reverse_map
 
     @property
-    def n_sample(self):
-        return self.shape[self.ax_sample] if self.ax_sample is not None else 1
+    def shape(self):
+        """Shape of array.  Unlike numpy ndarray, this is not assignable."""
+        return super().shape
+
+    def insert_axis(
+        self, axis: Union[Collection[int], int], new_name: str
+    ) -> CompatDict[int]:
+        """Create the constructor axes dict from this array, with new axis/axes"""
+        return self._ax_map.insert_axis(axis, new_name)
+
+    def remove_axis(self, axis: Union[Collection[int], int]) -> CompatDict[int]:
+        """Create the constructor axes dict from this array, without axis/axes"""
+        return self._ax_map.remove_axis(axis)
+
+    def __getattr__(self, name):
+        # TODO: replace with structural pattern matching on Oct 2025 (3.9 EOL)
+        parts = name.split("_", 1)
+        if parts[0] == "ax":
+            try:
+                return self.axes[name]
+            except KeyError:
+                raise AttributeError(f"AxesArray has no axis '{name}'")
+        if parts[0] == "n":
+            try:
+                ax_ids = self._ax_map.fwd_map["ax_" + parts[1]]
+            except KeyError:
+                raise AttributeError(f"AxesArray has no axis '{name}'")
+            shape = tuple(self.shape[ax_id] for ax_id in ax_ids)
+            if len(shape) == 1:
+                return shape[0]
+            return shape
+        raise AttributeError(f"'{type(self)}' object has no attribute '{name}'")
+
+    def __getitem__(self, key: Union[Indexer, Sequence[Indexer]], /):
+        if isinstance(key, tuple):
+            base_indexer = tuple(None if isinstance(k, str) else k for k in key)
+        else:
+            base_indexer = key
+        output = super().__getitem__(base_indexer)
+        if not isinstance(output, AxesArray):
+            return output  # return an element from the array
+        in_dim = self.shape
+        key, adv_inds = _standardize_indexer(self, key)
+        bcast_nd, bcast_start_ax = _determine_adv_broadcasting(key, adv_inds)
+        if adv_inds:
+            key = _replace_adv_indexers(key, adv_inds, bcast_start_ax, bcast_nd)
+        remove_axes, new_axes, adv_names = _apply_indexing(key, self._reverse_map)
+        new_axes = _rename_broadcast_axes(new_axes, adv_names)
+        new_map = _AxisMapping(
+            self._ax_map.remove_axis(remove_axes), len(in_dim) - len(remove_axes)
+        )
+        for insert_counter, (new_ax_ind, new_ax_name) in enumerate(new_axes):
+            new_map = _AxisMapping(
+                new_map.insert_axis(new_ax_ind, new_ax_name),
+                in_ndim=len(in_dim) - len(remove_axes) + (insert_counter + 1),
+            )
+        output._ax_map = new_map
+        return output
 
-    @property
-    def n_coord(self):
-        return self.shape[self.ax_coord] if self.ax_coord is not None else 1
+    def __array_finalize__(self, obj) -> None:
+        if obj is None:  # explicit construction via super().__new__()
+            return
+        # view from numpy array, called in constructor but also tests
+        if all(
+            (
+                not isinstance(obj, AxesArray),
+                self.shape == (),
+                not hasattr(self, "_ax_map"),
+            )
+        ):
+            self._ax_map = _AxisMapping({}, in_ndim=0)
+        # required by ravel() and view() used in numpy testing.  Also for zeros_like...
+        elif all(
+            (
+                isinstance(obj, AxesArray),
+                not hasattr(self, "_ax_map"),
+                self.shape == obj.shape,
+            )
+        ):
+            self._ax_map = _AxisMapping(obj.axes, obj.ndim)
+        # maybe add errors for incompatible views?
 
     def __array_ufunc__(
         self, ufunc, method, *inputs, out=None, **kwargs
@@ -87,27 +381,35 @@ def __array_ufunc__(
             return
         if ufunc.nout == 1:
             results = (results,)
-        results = tuple(
-            (AxesArray(np.asarray(result), self.__dict__) if output is None else output)
-            for result, output in zip(results, outputs)
-        )
+        if method == "reduce" and (
+            "keepdims" not in kwargs.keys() or kwargs["keepdims"] is False
+        ):
+            axes = None
+            if kwargs["axis"] is not None:
+                axes = self._ax_map.remove_axis(axis=kwargs["axis"])
+        else:
+            axes = self.axes
+        final_results = []
+        for result, output in zip(results, outputs):
+            if output is not None:
+                final_results.append(output)
+            elif axes is None:
+                final_results.append(result)
+            else:
+                final_results.append(AxesArray(np.asarray(result), axes))
+        results = tuple(final_results)
         return results[0] if len(results) == 1 else results
 
     def __array_function__(self, func, types, args, kwargs):
         if func not in HANDLED_FUNCTIONS:
-            arr = super(AxesArray, self).__array_function__(func, types, args, kwargs)
-            if isinstance(arr, np.ndarray):
-                return AxesArray(arr, axes=self.__dict__)
-            elif arr is not None:
-                return arr
-            return
+            return super(AxesArray, self).__array_function__(func, types, args, kwargs)
         if not all(issubclass(t, AxesArray) for t in types):
             return NotImplemented
         return HANDLED_FUNCTIONS[func](*args, **kwargs)
 
 
-def implements(numpy_function):
-    """Register an __array_function__ implementation for MyArray objects."""
+def _implements(numpy_function):
+    """Register an __array_function__ implementation for AxesArray objects."""
 
     def decorator(func):
         HANDLED_FUNCTIONS[numpy_function] = func
@@ -116,24 +418,381 @@ def decorator(func):
     return decorator
 
 
-@implements(np.concatenate)
+@_implements(np.ix_)
+def ix_(*args: AxesArray):
+    calc = np.ix_(*(np.asarray(arg) for arg in args))
+    ax_names = [list(arr.axes)[0] for arr in args]
+    axes = fwd_from_names(ax_names)
+    return tuple(AxesArray(arr, axes) for arr in calc)
+
+
+@_implements(np.concatenate)
 def concatenate(arrays, axis=0, out=None, dtype=None, casting="same_kind"):
     parents = [np.asarray(obj) for obj in arrays]
-    ax_list = [obj.__dict__ for obj in arrays if isinstance(obj, AxesArray)]
+    ax_list = [obj.axes for obj in arrays if isinstance(obj, AxesArray)]
     for ax1, ax2 in zip(ax_list[:-1], ax_list[1:]):
         if ax1 != ax2:
-            raise TypeError("Concatenating >1 AxesArray with incompatible axes")
+            raise ValueError("Concatenating >1 AxesArray with incompatible axes")
     result = np.concatenate(parents, axis, out=out, dtype=dtype, casting=casting)
     if isinstance(out, AxesArray):
-        out.__dict__ = ax_list[0]
+        out._ax_map = _AxisMapping(ax_list[0], in_ndim=result.ndim)
     return AxesArray(result, axes=ax_list[0])
 
 
+@_implements(np.reshape)
+def reshape(a: AxesArray, newshape: int | tuple[int], order="C"):
+    """Gives a new shape to an array without changing its data.
+
+    Args:
+        a: Array to be reshaped
+        newshape: int or tuple of ints
+            The new shape should be compatible with the original shape.  In
+            addition, the axis labels must make sense when the data is
+            translated to a new shape.  Currently, the only use case supported
+            is to flatten an outer product of two or more axes with the same
+            label and size.
+        order: Must be "C"
+    """
+    if order != "C":
+        raise ValueError("AxesArray only supports reshaping in 'C' order currently.")
+    out = np.reshape(np.asarray(a), newshape, order)  # handle any regular errors
+
+    new_axes = {}
+    if isinstance(newshape, int):
+        newshape = [newshape]
+    newshape = list(newshape)
+    explicit_new_size = np.multiply.reduce(np.array(newshape))
+    if explicit_new_size < 0:
+        replace_ind = newshape.index(-1)
+        newshape[replace_ind] = a.size // (-1 * explicit_new_size)
+
+    curr_base = 0
+    for curr_new in range(len(newshape)):
+        if curr_base >= a.ndim:
+            raise ValueError(
+                "Cannot reshape an AxesArray this way.  Adding a length-1 axis at"
+                f" dimension {curr_new} not understood."
+            )
+        base_name = a._ax_map.reverse_map[curr_base]
+        if a.shape[curr_base] == newshape[curr_new]:
+            compat_dict_append(new_axes, base_name, curr_new)
+            curr_base += 1
+        elif newshape[curr_new] == 1:
+            raise ValueError(
+                f"Cannot reshape an AxesArray this way.  Inserting a new axis at"
+                f" dimension {curr_new} of new shape is not supported"
+            )
+        else:  # outer product
+            remaining = newshape[curr_new]
+            while remaining > 1:
+                if a._ax_map.reverse_map[curr_base] != base_name:
+                    raise ValueError(
+                        "Cannot reshape an AxesArray this way.  It would combine"
+                        f" {base_name} with {a._ax_map.reverse_map[curr_base]}"
+                    )
+                remaining, error = divmod(remaining, a.shape[curr_base])
+                if error:
+                    raise ValueError(
+                        f"Cannot reshape an AxesArray this way.  Array dimension"
+                        f" {curr_base} has size {a.shape[curr_base]}, must divide into"
+                        f" newshape dimension {curr_new} with size"
+                        f" {newshape[curr_new]}."
+                    )
+                curr_base += 1
+
+            compat_dict_append(new_axes, base_name, curr_new)
+
+    return AxesArray(out, axes=new_axes)
+
+
+@_implements(np.transpose)
+def transpose(a: AxesArray, axes: Optional[Union[Tuple[int], List[int]]] = None):
+    """Returns an array with axes transposed.
+
+    Args:
+        a: input array
+        axes: As the numpy function
+    """
+    out = np.transpose(np.asarray(a), axes)
+    if axes is None:
+        axes = range(a.ndim)[::-1]
+    new_axes = {}
+    old_reverse = a._ax_map.reverse_map
+    for new_ind, old_ind in enumerate(axes):
+        compat_dict_append(new_axes, old_reverse[old_ind], new_ind)
+
+    return AxesArray(out, new_axes)
+
+
+@_implements(np.einsum)
+def einsum(
+    subscripts: str, *operands: AxesArray, out: Optional[NDArray] = None, **kwargs
+) -> AxesArray:
+    calc = np.einsum(
+        subscripts, *(np.asarray(arr) for arr in operands), out=out, **kwargs
+    )
+    try:
+        # explicit mode
+        lscripts, rscript = subscripts.split("->")
+    except ValueError:
+        # implicit mode
+        lscripts = subscripts
+        rscript = "".join(
+            sorted(c for c in set(subscripts) if subscripts.count(c) == 1 and c != ",")
+        )
+    # 0-dimensional case, may just be better to check type of "calc":
+    if rscript == "":
+        return calc
+
+    # assemble output reverse map
+    allscript_names = _label_einsum_scripts(lscripts, operands)
+    out_names = []
+
+    for char in rscript.replace("...", "."):
+        if char == ".":
+            for script_names in allscript_names:
+                out_names += script_names.get("...", [])
+        else:
+            ax_names = []
+            for script_names in allscript_names:
+                ax_names += script_names.get(char, [])
+            ax_name = "ax_" + _join_unique_names(ax_names)
+            out_names.append(ax_name)
+
+    out_axes = fwd_from_names(out_names)
+    if isinstance(out, AxesArray):
+        out._ax_map = _AxisMapping(out_axes, calc.ndim)
+    return AxesArray(calc, axes=out_axes)
+
+
+def _join_unique_names(l_of_s: List[str]) -> str:
+    ordered_uniques = dict.fromkeys(l_of_s).keys()
+    return "_".join(
+        ax_name[3:] if ax_name[:3] == "ax_" else ax_name for ax_name in ordered_uniques
+    )
+
+
+def _label_einsum_scripts(
+    lscripts: List[str], operands: tuple[AxesArray]
+) -> List[dict[str, str]]:
+    """Create a list of what axis name each script refers to in its operand."""
+    allscript_names: List[Dict[str, List[str]]] = []
+    for lscr, op in zip(lscripts.split(","), operands):
+        script_names: Dict[str, List[str]] = {}
+        allscript_names.append(script_names)
+        # handle script ellipses
+        try:
+            ell_ind = lscr.index("...")
+            ell_width = op.ndim - (len(lscr) - 3)
+            ell_expand = range(ell_ind, ell_ind + ell_width)
+            ell_names = [op._ax_map.reverse_map[ax_ind] for ax_ind in ell_expand]
+            script_names["..."] = ell_names
+        except ValueError:
+            ell_ind = len(lscr)
+            ell_width = 0
+        # handle script non-ellipsis chars
+        shift = 0
+        for ax_ind, char in enumerate(lscr):
+            if char == ".":
+                shift += 1
+                continue
+            if ax_ind < ell_ind:
+                scr_name = op._ax_map.reverse_map[ax_ind]
+            else:
+                scr_name = op._ax_map.reverse_map[ax_ind - 3 + ell_width]
+            compat_dict_append(script_names, char, [scr_name])
+    return allscript_names
+
+
+@_implements(np.linalg.solve)
+def linalg_solve(a: AxesArray, b: AxesArray) -> AxesArray:
+    result = np.linalg.solve(np.asarray(a), np.asarray(b))
+    a_rev = a._ax_map.reverse_map
+    a_names = [a_rev[k] for k in sorted(a_rev)]
+    contracted_axis_name = a_names[-1]
+    b_rev = b._ax_map.reverse_map
+    b_names = [b_rev[k] for k in sorted(b_rev)]
+    match_axes_list = a_names[:-1]
+    start = max(b.ndim - a.ndim, 0)
+    end = start + len(match_axes_list)
+    align = slice(start, end)
+    if match_axes_list != b_names[align]:
+        raise ValueError("Mismatch in operand axis names when aligning A and b")
+    all_names = (
+        b_names[: align.stop - 1] + [contracted_axis_name] + b_names[align.stop :]
+    )
+    axes = fwd_from_names(all_names)
+    return AxesArray(result, axes)
+
+
+@_implements(np.tensordot)
+def tensordot(
+    a: AxesArray, b: AxesArray, axes: Union[int, Sequence[Sequence[int]]] = 2
+) -> AxesArray:
+    sub = _tensordot_to_einsum(a.ndim, b.ndim, axes)
+    return einsum(sub, a, b)
+
+
+def _tensordot_to_einsum(
+    a_ndim: int, b_ndim: int, axes: Union[int, Sequence[Sequence[int]]]
+) -> str:
+    lc_ord = range(97, 123)
+    sub_a = "".join([chr(code) for code in lc_ord[:a_ndim]])
+    if isinstance(axes, int):
+        axes = [range(-axes, 0), range(0, axes)]
+    sub_b_li = [chr(code) for code in lc_ord[a_ndim : a_ndim + b_ndim]]
+    if np.array(axes).max() > 26:
+        raise ValueError("Too many axes")
+    for a_ind, b_ind in zip(*axes):
+        sub_b_li[b_ind] = sub_a[a_ind]
+    sub_b = "".join(sub_b_li)
+    sub = f"{sub_a},{sub_b}"
+    return sub
+
+
+def _standardize_indexer(
+    arr: np.ndarray, key: Indexer | Sequence[Indexer]
+) -> tuple[Sequence[StandardIndexer], tuple[KeyIndex, ...]]:
+    """Convert any legal numpy indexer to a "standard" form.
+
+    Standard form involves creating an equivalent indexer that is a tuple with
+    one element per index of the original axis.  All advanced indexer elements
+    are converted to numpy arrays, and boolean arrays are converted to
+    integer arrays with obj.nonzero().
+
+    Returns:
+        A tuple of the normalized indexer as well as the indexes of
+        advanced indexers
+    """
+    if isinstance(key, tuple):
+        key = list(key)
+    else:
+        key = [key]
+
+    if not any(ax_key is Ellipsis for ax_key in key):
+        key = [*key, Ellipsis]
+
+    new_key: List[Indexer] = []
+    for ax_key in key:
+        if not isinstance(ax_key, get_args(BasicIndexer)):
+            ax_key = np.array(ax_key)
+            if ax_key.dtype == np.dtype(np.bool_):
+                new_key += ax_key.nonzero()
+                continue
+        new_key.append(ax_key)
+
+    new_key = _expand_indexer_ellipsis(new_key, arr.ndim)
+    # Can't identify position of advanced indexers before expanding ellipses
+    adv_inds: List[KeyIndex] = []
+    for key_ind, ax_key in enumerate(new_key):
+        if isinstance(ax_key, np.ndarray):
+            adv_inds.append(KeyIndex(key_ind))
+
+    return new_key, tuple(adv_inds)
+
+
+def _expand_indexer_ellipsis(key: List[Indexer], ndim: int) -> List[Indexer]:
+    """Replace ellipsis in indexers with the appropriate amount of slice(None)"""
+    # [...].index errors if list contains numpy array
+    ellind = [ind for ind, val in enumerate(key) if val is ...][0]
+    n_new_dims = sum(ax_key is None or isinstance(ax_key, str) for ax_key in key)
+    n_ellipsis_dims = ndim - (len(key) - n_new_dims - 1)
+    new_key = key[:ellind] + key[ellind + 1 :]
+    new_key = new_key[:ellind] + (n_ellipsis_dims * [slice(None)]) + new_key[ellind:]
+    return new_key
+
+
+def _determine_adv_broadcasting(
+    key: Sequence[StandardIndexer], adv_inds: Sequence[OldIndex]
+) -> tuple[int, Optional[KeyIndex]]:
+    """Calculate the shape and location for the result of advanced indexing."""
+    adjacent = all(i + 1 == j for i, j in zip(adv_inds[:-1], adv_inds[1:]))
+    adv_indexers = [np.array(key[i]) for i in adv_inds]
+    bcast_nd = np.broadcast(*adv_indexers).nd
+    bcast_start_axis = 0 if not adjacent else min(adv_inds) if adv_inds else None
+    return bcast_nd, KeyIndex(bcast_start_axis)
+
+
+def _rename_broadcast_axes(
+    new_axes: List[tuple[int, None | str | Literal[_Sentinels.ADV_NAME]]],
+    adv_names: List[str],
+) -> List[tuple[int, str]]:
+    """Normalize sentinel and NoneType names"""
+
+    def _calc_bcast_name(*names: str) -> str:
+        if not names:
+            return ""
+        if all(a == b for a, b in zip(names[1:], names[:-1])):
+            return names[0]
+        names = [name[3:] for name in dict.fromkeys(names)]  # ordered deduplication
+        return "ax_" + "_".join(names)
+
+    bcast_name = _calc_bcast_name(*adv_names)
+    renamed_axes = []
+    for ax_ind, ax_name in new_axes:
+        if ax_name is None:
+            renamed_axes.append((ax_ind, "ax_unk"))
+        elif ax_name is _Sentinels.ADV_NAME:
+            renamed_axes.append((ax_ind, bcast_name))
+        else:
+            renamed_axes.append((ax_ind, "ax_" + ax_name))
+    return renamed_axes
+
+
+def _replace_adv_indexers(
+    key: Sequence[StandardIndexer],
+    adv_inds: List[int],
+    bcast_start_ax: int,
+    bcast_nd: int,
+) -> tuple[
+    Union[None, str, int, Literal[_Sentinels.ADV_NAME], Literal[_Sentinels.ADV_REMOVE]],
+    ...,
+]:
+    for adv_ind in adv_inds:
+        key[adv_ind] = _Sentinels.ADV_REMOVE
+    key = key[:bcast_start_ax] + bcast_nd * [_Sentinels.ADV_NAME] + key[bcast_start_ax:]
+    return key
+
+
+def _apply_indexing(
+    key: tuple[StandardIndexer], reverse_map: Dict[int, str]
+) -> tuple[
+    List[int], List[tuple[int, None | str | Literal[_Sentinels.ADV_NAME]]], List[str]
+]:
+    """Determine where axes should be removed and added
+
+    Only considers the basic indexers in key.  Numpy arrays are treated as
+    slices, in that they don't affect the final dimensions of the output
+    """
+    remove_axes = []
+    new_axes = []
+    adv_names = []
+    deleted_to_left = 0
+    added_to_left = 0
+    for key_ind, indexer in enumerate(key):
+        if isinstance(indexer, int) or indexer is _Sentinels.ADV_REMOVE:
+            orig_arr_axis = key_ind - added_to_left
+            if indexer is _Sentinels.ADV_REMOVE:
+                adv_names.append(reverse_map[orig_arr_axis])
+            remove_axes.append(orig_arr_axis)
+            deleted_to_left += 1
+        elif (
+            indexer is None
+            or indexer is _Sentinels.ADV_NAME
+            or isinstance(indexer, str)
+        ):
+            new_arr_axis = key_ind - deleted_to_left
+            new_axes.append((new_arr_axis, indexer))
+            added_to_left += 1
+    return remove_axes, new_axes, adv_names
+
+
 def comprehend_axes(x):
     axes = {}
     axes["ax_coord"] = len(x.shape) - 1
     axes["ax_time"] = len(x.shape) - 2
-    axes["ax_spatial"] = list(range(len(x.shape) - 2))
+    if x.ndim > 2:
+        axes["ax_spatial"] = list(range(len(x.shape) - 2))
     return axes
 
 
@@ -155,13 +814,15 @@ def concat_sample_axis(x_list: List[AxesArray]):
     """Concatenate all trajectories and axes used to create samples."""
     new_arrs = []
     for x in x_list:
-        sample_axes = (
-            x.ax_spatial
-            + ([x.ax_time] if x.ax_time is not None else [])
-            + ([x.ax_sample] if x.ax_sample is not None else [])
-        )
+        sample_ax_names = ("ax_spatial", "ax_time", "ax_sample")
+        sample_ax_inds = []
+        for name in sample_ax_names:
+            ax_inds = getattr(x, name, [])
+            if isinstance(ax_inds, int):
+                ax_inds = [ax_inds]
+            sample_ax_inds += ax_inds
         new_axes = {"ax_sample": 0, "ax_coord": 1}
-        n_samples = np.prod([x.shape[ax] for ax in sample_axes])
+        n_samples = np.prod([x.shape[ax] for ax in sample_ax_inds])
         arr = AxesArray(x.reshape((n_samples, x.shape[x.ax_coord])), new_axes)
         new_arrs.append(arr)
     return np.concatenate(new_arrs, axis=new_arrs[0].ax_sample)
@@ -176,3 +837,29 @@ def wrap_axes(axes: dict, obj):
         except KeyError:
             pass
     return obj
+
+
+def compat_dict_append(
+    compat_dict: CompatDict[T],
+    key: str,
+    item_or_list: ItemOrList[T],
+) -> None:
+    """Add an element or list of elements to a dictionary, preserving old values"""
+    try:
+        prev_val = compat_dict[key]
+    except KeyError:
+        compat_dict[key] = item_or_list
+        return
+    if not isinstance(item_or_list, list):
+        item_or_list = [item_or_list]
+    if not isinstance(prev_val, list):
+        prev_val = [prev_val]
+    compat_dict[key] = prev_val + item_or_list
+
+
+def fwd_from_names(names: List[str]) -> CompatDict[int]:
+    """Create mapping of name: axis or name: [ax_1, ax_2, ...]"""
+    fwd_map: Dict[str, Sequence[int]] = {}
+    for ax_ind, name in enumerate(names):
+        compat_dict_append(fwd_map, name, [ax_ind])
+    return fwd_map
diff --git a/test/test_feature_library.py b/test/test_feature_library.py
index bd5a05d2c..3b2f032cb 100644
--- a/test/test_feature_library.py
+++ b/test/test_feature_library.py
@@ -235,6 +235,7 @@ def test_sindypi_library_bad_params(params):
         pytest.lazy_fixture("ode_library"),
         pytest.lazy_fixture("sindypi_library"),
     ],
+    ids=type,
 )
 def test_fit_transform(data_lorenz, library):
     x, t = data_lorenz
diff --git a/test/test_optimizers.py b/test/test_optimizers.py
index 171825c4b..f59b7dd47 100644
--- a/test/test_optimizers.py
+++ b/test/test_optimizers.py
@@ -589,7 +589,7 @@ def test_specific_bad_parameters(error, optimizer, params, data_lorenz):
 def test_bad_optimizers(data_derivative_1d):
     x, x_dot = data_derivative_1d
     x = x.reshape(-1, 1)
-
+    x_dot = x_dot.reshape(-1, 1)
     with pytest.raises(InvalidParameterError):
         # Error: optimizer does not have a callable fit method
         opt = WrappedOptimizer(DummyEmptyModel())
diff --git a/test/test_optimizers_complexity.py b/test/test_optimizers_complexity.py
index 8a6486d83..bcdbe522c 100644
--- a/test/test_optimizers_complexity.py
+++ b/test/test_optimizers_complexity.py
@@ -12,6 +12,7 @@
 from pysindy.optimizers import WrappedOptimizer
 
 
+@pytest.mark.skip
 @pytest.mark.parametrize(
     "opt_cls, reg_name",
     [[Lasso, "alpha"], [Ridge, "alpha"], [STLSQ, "threshold"], [SR3, "threshold"]],
diff --git a/test/utils/test_axes.py b/test/utils/test_axes.py
index e5d9a8385..b26a73890 100644
--- a/test/utils/test_axes.py
+++ b/test/utils/test_axes.py
@@ -1,10 +1,18 @@
 import numpy as np
 import pytest
 from numpy.testing import assert_
+from numpy.testing import assert_array_equal
 from numpy.testing import assert_equal
 from numpy.testing import assert_raises
 
 from pysindy import AxesArray
+from pysindy.utils import axes
+from pysindy.utils.axes import _AxisMapping
+from pysindy.utils.axes import AxesWarning
+
+
+def test_axesarray_create():
+    AxesArray(np.array(1), {})
 
 
 def test_concat_out():
@@ -14,11 +22,25 @@ def test_concat_out():
     assert_equal(result, arr_out)
 
 
+def test_bad_concat():
+    arr = AxesArray(np.arange(3).reshape(1, 3), {"ax_a": 0, "ax_b": 1})
+    arr2 = AxesArray(np.arange(3).reshape(1, 3), {"ax_b": 0, "ax_c": 1})
+    with pytest.raises(ValueError):
+        np.concatenate((arr, arr2), axis=0)
+
+
 def test_reduce_mean_noinf_recursion():
-    arr = AxesArray(np.array([[1]]), {})
+    arr = AxesArray(np.array([[1]]), {"ax_a": [0, 1]})
     np.mean(arr, axis=0)
 
 
+def test_repr():
+    a = AxesArray(np.arange(5.0), {"ax_time": 0})
+    result = a.__repr__()
+    expected = "AxesArray([0., 1., 2., 3., 4.])"
+    assert result == expected
+
+
 def test_ufunc_override():
     # This is largely a clone of test_ufunc_override_with_super() from
     # numpy/core/tests/test_umath.py
@@ -32,31 +54,31 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
 
     d = np.arange(5.0)
     # 1 input, 1 output
-    a = AxesArray(d, {})
+    a = AxesArray(d, {"ax_time": 0})
     b = np.sin(a)
     check = np.sin(d)
     assert_(np.all(check == b))
     b = np.sin(d, out=(a,))
     assert_(np.all(check == b))
     assert_(b is a)
-    a = AxesArray(np.arange(5.0), {})
+    a = AxesArray(np.arange(5.0), {"ax_time": 0})
     b = np.sin(a, out=a)
     assert_(np.all(check == b))
 
     # 1 input, 2 outputs
-    a = AxesArray(np.arange(5.0), {})
+    a = AxesArray(np.arange(5.0), {"ax_time": 0})
     b1, b2 = np.modf(a)
     b1, b2 = np.modf(d, out=(None, a))
     assert_(b2 is a)
-    a = AxesArray(np.arange(5.0), {})
-    b = AxesArray(np.arange(5.0), {})
+    a = AxesArray(np.arange(5.0), {"ax_time": 0})
+    b = AxesArray(np.arange(5.0), {"ax_time": 0})
     c1, c2 = np.modf(a, out=(a, b))
     assert_(c1 is a)
     assert_(c2 is b)
 
     # 2 input, 1 output
-    a = AxesArray(np.arange(5.0), {})
-    b = AxesArray(np.arange(5.0), {})
+    a = AxesArray(np.arange(5.0), {"ax_time": 0})
+    b = AxesArray(np.arange(5.0), {"ax_time": 0})
     c = np.add(a, b, out=a)
     assert_(c is a)
     # some tests with a non-ndarray subclass
@@ -65,13 +87,13 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
     assert_(a.__array_ufunc__(np.add, "__call__", a, b) is NotImplemented)
     assert_(b.__array_ufunc__(np.add, "__call__", a, b) is NotImplemented)
     assert_raises(TypeError, np.add, a, b)
-    a = AxesArray(a, {})
+    a = AxesArray(a, {"ax_time": 0})
     assert_(a.__array_ufunc__(np.add, "__call__", a, b) is NotImplemented)
     assert_(b.__array_ufunc__(np.add, "__call__", a, b) == "A!")
     assert_(np.add(a, b) == "A!")
     # regression check for gh-9102 -- tests ufunc.reduce implicitly.
     d = np.array([[1, 2, 3], [1, 2, 3]])
-    a = AxesArray(d, {})
+    a = AxesArray(d, {"ax_time": [0, 1]})
     c = a.any()
     check = d.any()
     assert_equal(c, check)
@@ -116,11 +138,11 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
     a = d.copy().view(AxesArray)
     np.add.at(check, ([0, 1], [0, 2]), 1.0)
     np.add.at(a, ([0, 1], [0, 2]), 1.0)
-    assert_equal(a, check)
+    assert_equal(np.asarray(a), np.asarray(check))  # modified
     b = np.array(1.0).view(AxesArray)
     a = d.copy().view(AxesArray)
     np.add.at(a, ([0, 1], [0, 2]), b)
-    assert_equal(a, check)
+    assert_equal(np.asarray(a), np.asarray(check))  # modified
 
 
 def test_n_elements():
@@ -129,18 +151,526 @@ def test_n_elements():
     assert arr.n_spatial == (1, 2)
     assert arr.n_time == 3
     assert arr.n_coord == 4
-    assert arr.n_sample == 1
 
     arr2 = np.concatenate((arr, arr), axis=arr.ax_time)
     assert arr2.n_spatial == (1, 2)
     assert arr2.n_time == 6
     assert arr2.n_coord == 4
-    assert arr2.n_sample == 1
-
-    arr3 = arr[..., :2, 0]
-    assert arr3.n_spatial == (1, 2)
-    assert arr3.n_time == 2
-    # No way to intercept slicing and remove ax_coord
-    with pytest.raises(IndexError):
-        assert arr3.n_coord == 1
-    assert arr3.n_sample == 1
+
+
+def test_reshape_outer_product():
+    arr = AxesArray(np.arange(4).reshape((2, 2)), {"ax_a": [0, 1]})
+    merge = np.reshape(arr, (4,))
+    assert merge.axes == {"ax_a": 0}
+
+
+def test_reshape_bad_divmod():
+    arr = AxesArray(np.arange(12).reshape((2, 3, 2)), {"ax_a": [0, 1], "ax_b": 2})
+    with pytest.raises(
+        ValueError, match="Cannot reshape an AxesArray this way.  Array dimension"
+    ):
+        np.reshape(arr, (4, 3))
+
+
+def test_reshape_fill_outer_product():
+    arr = AxesArray(np.arange(4).reshape((2, 2)), {"ax_a": [0, 1]})
+    merge = np.reshape(arr, (-1,))
+    assert merge.axes == {"ax_a": 0}
+
+
+def test_reshape_fill_regular():
+    arr = AxesArray(np.arange(8).reshape((2, 2, 2)), {"ax_a": [0, 1], "ax_b": 2})
+    merge = np.reshape(arr, (4, -1))
+    assert merge.axes == {"ax_a": 0, "ax_b": 1}
+
+
+def test_illegal_reshape():
+    arr = AxesArray(np.arange(4).reshape((2, 2)), {"ax_a": [0, 1]})
+    # melding across axes
+    with pytest.raises(ValueError, match="Cannot reshape an AxesArray"):
+        np.reshape(arr, (4, 1))
+
+    # Add a hidden 1 in the middle!  maybe a matching 1
+
+    # different name outer product
+    arr = AxesArray(np.arange(4).reshape((2, 2)), {"ax_a": 0, "ax_b": 1})
+    with pytest.raises(ValueError, match="Cannot reshape an AxesArray"):
+        np.reshape(arr, (4,))
+    # newaxes
+    with pytest.raises(ValueError, match="Cannot reshape an AxesArray"):
+        np.reshape(arr, (2, 1, 2))
+
+
+def test_warn_toofew_axes():
+    axes = {"ax_time": 0, "ax_coord": 1}
+    with pytest.warns(AxesWarning):
+        AxesArray(np.ones(8).reshape((2, 2, 2)), axes)
+
+
+def test_toomany_axes():
+    axes = {"ax_time": 0, "ax_coord": 2}
+    with pytest.raises(ValueError):
+        AxesArray(np.ones(4).reshape((2, 2)), axes)
+
+
+def test_conflicting_axes_defn():
+    axes = {"ax_time": 0, "ax_coord": 0}
+    with pytest.raises(ValueError):
+        AxesArray(np.ones(4), axes)
+
+
+def test_missing_axis_errors():
+    axes = {"ax_time": 0}
+    arr = AxesArray(np.arange(3), axes)
+    with pytest.raises(AttributeError):
+        arr.ax_spatial
+    with pytest.raises(AttributeError):
+        arr.n_spatial
+
+
+def test_simple_slice():
+    arr = AxesArray(np.ones(2), {"ax_coord": 0})
+    assert_array_equal(arr[:], arr)
+    assert_array_equal(arr[slice(None)], arr)
+    assert arr[0] == 1
+
+
+# @pytest.mark.skip  # TODO: make this pass
+def test_0d_indexer():
+    arr = AxesArray(np.ones(2), {"ax_coord": 0})
+    arr_out = arr[1, ...]
+    assert arr_out.ndim == 0
+    assert arr_out.axes == {}
+    assert arr_out[()] == 1
+
+
+def test_basic_indexing_modifies_axes():
+    axes = {"ax_time": 0, "ax_coord": 1}
+    arr = AxesArray(np.ones(4).reshape((2, 2)), axes)
+    slim = arr[1, :, None]
+    with pytest.raises(AttributeError):
+        slim.ax_time
+    assert slim.ax_unk == 1
+    assert slim.ax_coord == 0
+    reverse_slim = arr[None, :, 1]
+    with pytest.raises(AttributeError):
+        reverse_slim.ax_coord
+    assert reverse_slim.ax_unk == 0
+    assert reverse_slim.ax_time == 1
+    almost_new = arr[None, None, 1, :, None, None]
+    with pytest.raises(AttributeError):
+        almost_new.ax_time
+    assert almost_new.ax_coord == 2
+    assert set(almost_new.ax_unk) == {0, 1, 3, 4}
+
+
+def test_insert_named_axis():
+    arr = AxesArray(np.ones(1), axes={"ax_time": 0})
+    expanded = arr["time", :]
+    result = expanded.axes
+    expected = {"ax_time": [0, 1]}
+    assert result == expected
+
+
+def test_adv_indexing_modifies_axes():
+    axes = {"ax_time": 0, "ax_coord": 1}
+    arr = AxesArray(np.arange(4).reshape((2, 2)), axes)
+    flat = arr[[0, 1], [0, 1]]
+    same = arr[[[0], [1]], [0, 1]]
+    tpose = arr[[0, 1], [[0], [1]]]
+    assert flat.shape == (2,)
+    np.testing.assert_array_equal(np.asarray(flat), np.array([0, 3]))
+
+    assert flat.ax_time_coord == 0
+    with pytest.raises(AttributeError):
+        flat.ax_coord
+    with pytest.raises(AttributeError):
+        flat.ax_time
+
+    assert same.shape == arr.shape
+    np.testing.assert_equal(np.asarray(same), np.asarray(arr))
+    assert same.ax_time_coord == [0, 1]
+    with pytest.raises(AttributeError):
+        same.ax_coord
+
+    assert tpose.shape == arr.shape
+    np.testing.assert_equal(np.asarray(tpose), np.asarray(arr.T))
+    assert tpose.ax_time_coord == [0, 1]
+    with pytest.raises(AttributeError):
+        tpose.ax_coord
+
+
+def test_adv_indexing_adds_axes():
+    axes = {"ax_time": 0, "ax_coord": 1}
+    arr = AxesArray(np.arange(4).reshape((2, 2)), axes)
+    fat = arr[[[0, 1], [0, 1]]]
+    assert fat.shape == (2, 2, 2)
+    assert fat.ax_time == [0, 1]
+    assert fat.ax_coord == 2
+
+
+def test_standardize_basic_indexer():
+    arr = np.arange(6).reshape(2, 3)
+    result_indexer, result_fancy = axes._standardize_indexer(arr, Ellipsis)
+    assert result_indexer == [slice(None), slice(None)]
+    assert result_fancy == ()
+
+    result_indexer, result_fancy = axes._standardize_indexer(
+        arr, (np.newaxis, 1, 1, Ellipsis)
+    )
+    assert result_indexer == [None, 1, 1]
+    assert result_fancy == ()
+
+
+def test_standardize_advanced_indexer():
+    arr = np.arange(6).reshape(2, 3)
+    result_indexer, result_fancy = axes._standardize_indexer(arr, [1])
+    assert result_indexer == [np.ones(1), slice(None)]
+    assert result_fancy == (0,)
+
+    result_indexer, result_fancy = axes._standardize_indexer(
+        arr, (np.newaxis, [1], 1, Ellipsis)
+    )
+    assert result_indexer == [None, np.ones(1), 1]
+    assert result_fancy == (1,)
+
+
+def test_standardize_bool_indexer():
+    arr = np.ones((1, 2))
+    result, result_adv = axes._standardize_indexer(arr, [[True, True]])
+    assert_equal(result, [[0, 0], [0, 1]])
+    assert result_adv == (0, 1)
+
+
+def test_reduce_AxisMapping():
+    ax_map = _AxisMapping(
+        {"ax_a": [0, 1], "ax_b": 2, "ax_c": 3, "ax_d": 4, "ax_e": [5, 6]},
+        7,
+    )
+    result = ax_map.remove_axis(3)
+    expected = {"ax_a": [0, 1], "ax_b": 2, "ax_d": 3, "ax_e": [4, 5]}
+    assert result == expected
+    result = ax_map.remove_axis(-4)
+    assert result == expected
+
+
+def test_reduce_all_AxisMapping():
+    ax_map = _AxisMapping({"ax_a": [0, 1], "ax_b": 2}, 3)
+    result = ax_map.remove_axis()
+    expected = {}
+    assert result == expected
+
+
+def test_reduce_multiple_AxisMapping():
+    ax_map = _AxisMapping(
+        {
+            "ax_a": [0, 1],
+            "ax_b": 2,
+            "ax_c": 3,
+            "ax_d": 4,
+            "ax_e": [5, 6],
+        },
+        7,
+    )
+    result = ax_map.remove_axis([3, 4])
+    expected = {
+        "ax_a": [0, 1],
+        "ax_b": 2,
+        "ax_e": [3, 4],
+    }
+    assert result == expected
+
+
+def test_reduce_twisted_AxisMapping():
+    ax_map = _AxisMapping(
+        {
+            "ax_a": [0, 6],
+            "ax_b": 2,
+            "ax_c": 3,
+            "ax_d": 4,
+            "ax_e": [1, 5],
+        },
+        7,
+    )
+    result = ax_map.remove_axis([3, 4])
+    expected = {
+        "ax_a": [0, 4],
+        "ax_b": 2,
+        "ax_e": [1, 3],
+    }
+    assert result == expected
+
+
+def test_reduce_misordered_AxisMapping():
+    ax_map = _AxisMapping({"ax_a": [0, 1], "ax_b": 2, "ax_c": 3}, 4)
+    result = ax_map.remove_axis([2, 1])
+    expected = {"ax_a": 0, "ax_c": 1}
+    assert result == expected
+
+
+def test_insert_AxisMapping():
+    ax_map = _AxisMapping(
+        {
+            "ax_a": [0, 1],
+            "ax_b": 2,
+            "ax_c": 3,
+            "ax_d": [4, 5],
+        },
+        6,
+    )
+    result = ax_map.insert_axis(3, "ax_unk")
+    expected = {
+        "ax_a": [0, 1],
+        "ax_b": 2,
+        "ax_unk": 3,
+        "ax_c": 4,
+        "ax_d": [5, 6],
+    }
+    assert result == expected
+
+
+def test_insert_existing_AxisMapping():
+    ax_map = _AxisMapping(
+        {
+            "ax_a": [0, 1],
+            "ax_b": 2,
+            "ax_c": 3,
+            "ax_d": [4, 5],
+        },
+        6,
+    )
+    result = ax_map.insert_axis(3, "ax_b")
+    expected = {
+        "ax_a": [0, 1],
+        "ax_b": [2, 3],
+        "ax_c": 4,
+        "ax_d": [5, 6],
+    }
+    assert result == expected
+
+
+def test_insert_multiple_AxisMapping():
+    ax_map = _AxisMapping(
+        {
+            "ax_a": [0, 1],
+            "ax_b": 2,
+            "ax_c": 3,
+            "ax_d": [4, 5],
+        },
+        6,
+    )
+    result = ax_map.insert_axis([1, 4], new_name="ax_unk")
+    expected = {
+        "ax_a": [0, 2],
+        "ax_unk": [1, 4],
+        "ax_b": 3,
+        "ax_c": 5,
+        "ax_d": [6, 7],
+    }
+    assert result == expected
+
+
+def test_insert_misordered_AxisMapping():
+    ax_map = _AxisMapping(
+        {
+            "ax_a": [0, 1],
+            "ax_b": 2,
+            "ax_c": 3,
+            "ax_d": [4, 5],
+        },
+        6,
+    )
+    result = ax_map.insert_axis([4, 1], new_name="ax_unk")
+    expected = {
+        "ax_a": [0, 2],
+        "ax_unk": [1, 4],
+        "ax_b": 3,
+        "ax_c": 5,
+        "ax_d": [6, 7],
+    }
+    assert result == expected
+
+
+def test_determine_adv_broadcasting():
+    indexers = (1, np.ones(1), np.ones((4, 1)), np.ones(3))
+    res_nd, res_start = axes._determine_adv_broadcasting(indexers, [1, 2, 3])
+    assert res_nd == 2
+    assert res_start == 1
+
+    indexers = (None, np.ones(1), 2, np.ones(3))
+    res_nd, res_start = axes._determine_adv_broadcasting(indexers, [1, 3])
+    assert res_nd == 1
+    assert res_start == 0
+
+    res_nd, res_start = axes._determine_adv_broadcasting(indexers, [])
+    assert res_nd == 0
+    assert res_start is None
+
+
+def test_replace_ellipsis():
+    key = [..., 0]
+    result = axes._expand_indexer_ellipsis(key, 2)
+    expected = [slice(None), 0]
+    assert result == expected
+
+
+def test_strip_ellipsis():
+    key = [1, ...]
+    result = axes._expand_indexer_ellipsis(key, 1)
+    expected = [1]
+    assert result == expected
+
+    key = [..., 1]
+    result = axes._expand_indexer_ellipsis(key, 1)
+    expected = [1]
+    assert result == expected
+
+
+def test_transpose():
+    axes = {"ax_a": 0, "ax_b": [1, 2]}
+    arr = AxesArray(np.arange(8).reshape(2, 2, 2), axes)
+    tp = np.transpose(arr, [2, 0, 1])
+    result = tp.axes
+    expected = {"ax_a": 1, "ax_b": [0, 2]}
+    assert result == expected
+    assert_array_equal(tp, np.transpose(np.asarray(arr), [2, 0, 1]))
+    arr = arr[..., 0]
+    tp = arr.T
+    expected = {"ax_a": 1, "ax_b": 0}
+    assert_array_equal(tp, np.asarray(arr).T)
+
+
+def test_linalg_solve_align_left():
+    axesA = {"ax_prob": 0, "ax_sample": 1, "ax_coord": 2}
+    arrA = AxesArray(np.arange(8).reshape(2, 2, 2), axesA)
+    axesb = {"ax_prob": 0, "ax_sample": 1}
+    arrb = AxesArray(np.arange(4).reshape(2, 2), axesb)
+    result = np.linalg.solve(arrA, arrb)
+    expected_axes = {"ax_prob": 0, "ax_coord": 1}
+    assert result.axes == expected_axes
+    super_result = np.linalg.solve(np.asarray(arrA), np.asarray(arrb))
+    assert_array_equal(result, super_result)
+
+
+def test_linalg_solve_align_right():
+    axesA = {"ax_sample": 0, "ax_feature": 1}
+    arrA = AxesArray(np.arange(4).reshape(2, 2), axesA)
+    axesb = {"ax_sample": 0, "ax_target": 1}
+    arrb = AxesArray(np.arange(4).reshape(2, 2), axesb)
+    result = np.linalg.solve(arrA, arrb)
+    expected_axes = {"ax_feature": 0, "ax_target": 1}
+    assert result.axes == expected_axes
+    super_result = np.linalg.solve(np.asarray(arrA), np.asarray(arrb))
+    assert_array_equal(result, super_result)
+
+
+def test_linalg_solve_align_right_xl():
+    axesA = {"ax_sample": 0, "ax_feature": 1}
+    arrA = AxesArray(np.arange(4).reshape(2, 2), axesA)
+    axesb = {"ax_prob": 0, "ax_sample": 1, "ax_target": 2}
+    arrb = AxesArray(np.arange(8).reshape(2, 2, 2), axesb)
+    result = np.linalg.solve(arrA, arrb)
+    expected_axes = {"ax_prob": 0, "ax_feature": 1, "ax_target": 2}
+    assert result.axes == expected_axes
+    super_result = np.linalg.solve(np.asarray(arrA), np.asarray(arrb))
+    assert_array_equal(result, super_result)
+
+
+def test_linalg_solve_incompatible_left():
+    axesA = {"ax_prob": 0, "ax_sample": 1, "ax_coord": 2}
+    arrA = AxesArray(np.arange(8).reshape(2, 2, 2), axesA)
+    axesb = {"ax_foo": 0, "ax_sample": 1}
+    arrb = AxesArray(np.arange(4).reshape(2, 2), axesb)
+    with pytest.raises(ValueError, match="Mismatch in operand axis names"):
+        np.linalg.solve(arrA, arrb)
+
+
+def test_ts_to_einsum_int_axes():
+    a_str, b_str = axes._tensordot_to_einsum(3, 3, 2).split(",")
+    # expecting 'abc,bcf
+    assert a_str[0] not in b_str
+    assert b_str[-1] not in a_str
+    assert a_str[1:] == b_str[:-1]
+
+
+def test_ts_to_einsum_list_axes():
+    a_str, b_str = axes._tensordot_to_einsum(3, 3, [[1], [2]]).split(",")
+    # expecting 'abcd,efbh
+    assert a_str[1] == b_str[2]
+    assert a_str[0] not in b_str
+    assert a_str[2] not in b_str
+    assert b_str[0] not in a_str
+    assert b_str[1] not in a_str
+
+
+def test_tensordot_int_axes():
+    axes_a = {"ax_a": 0, "ax_b": [1, 2]}
+    axes_b = {"ax_b": [0, 1], "ax_c": 2}
+    arr = np.arange(8).reshape((2, 2, 2))
+    arr_a = AxesArray(arr, axes_a)
+    arr_b = AxesArray(arr, axes_b)
+    super_result = np.tensordot(arr, arr, 2)
+    result = np.tensordot(arr_a, arr_b, 2)
+    expected_axes = {"ax_a": 0, "ax_c": 1}
+    assert result.axes == expected_axes
+    assert_array_equal(result, super_result)
+
+
+def test_tensordot_list_axes():
+    axes_a = {"ax_a": 0, "ax_b": [1, 2]}
+    axes_b = {"ax_c": [0, 1], "ax_b": 2}
+    arr = np.arange(8).reshape((2, 2, 2))
+    arr_a = AxesArray(arr, axes_a)
+    arr_b = AxesArray(arr, axes_b)
+    super_result = np.tensordot(arr, arr, [[1], [2]])
+    result = np.tensordot(arr_a, arr_b, [[1], [2]])
+    expected_axes = {"ax_a": 0, "ax_b": 1, "ax_c": [2, 3]}
+    assert result.axes == expected_axes
+    assert_array_equal(result, super_result)
+
+
+@pytest.mark.skip
+def test_einsum_implicit():
+    ...
+
+
+@pytest.mark.skip
+def test_einsum_trace():
+    ...
+
+
+@pytest.mark.skip
+def test_einsum_diag():
+    ...
+
+
+@pytest.mark.skip
+def test_einsum_1dsum():
+    ...
+
+
+@pytest.mark.skip
+def test_einsum_alldsum():
+    ...
+
+
+@pytest.mark.skip
+def test_einsum_contraction():
+    ...
+
+
+@pytest.mark.skip
+def test_einsum_explicit_ellipsis():
+    ...
+
+
+def test_einsum_scalar():
+    arr = AxesArray(np.ones(1), {"ax_a": 0})
+    expected = 1
+    result = np.einsum("i,i", arr, arr)
+    assert result == expected
+
+
+@pytest.mark.skip
+def test_einsum_mixed():
+    ...