diff --git a/burnman/classes/anisotropy.py b/burnman/classes/anisotropy.py
index 2d945d9f..2f8cf70e 100644
--- a/burnman/classes/anisotropy.py
+++ b/burnman/classes/anisotropy.py
@@ -167,7 +167,7 @@ def isentropic_universal_elastic_anisotropy(self):
     def isentropic_isotropic_poisson_ratio(self):
         """
         :returns: The isotropic Poisson ratio (mu) [unitless].
-            A metric of the laterial response to loading.
+            A metric of the lateral response to loading.
         :rtype: float
         """
         return (
@@ -187,14 +187,11 @@ def christoffel_tensor(self, propagation_direction):
         :rtype: float
         """
         propagation_direction = unit_normalize(propagation_direction)
-        Tik = np.tensordot(
-            np.tensordot(
-                self.full_isentropic_stiffness_tensor,
-                propagation_direction,
-                axes=([1], [0]),
-            ),
+        Tik = np.einsum(
+            "ijkl, ...j, ...l",
+            self.full_isentropic_stiffness_tensor,
+            propagation_direction,
             propagation_direction,
-            axes=([2], [0]),
         )
         return Tik
 
@@ -205,8 +202,12 @@ def isentropic_linear_compressibility(self, direction):
         :rtype: float
         """
         direction = unit_normalize(direction)
-        Sijkk = np.einsum("ijkk", self.full_isentropic_compliance_tensor)
-        beta = Sijkk.dot(direction).dot(direction)
+        beta = np.einsum(
+            "ijkk, ...i, ...j",
+            self.full_isentropic_compliance_tensor,
+            direction,
+            direction,
+        )
         return beta
 
     def isentropic_youngs_modulus(self, direction):
@@ -216,8 +217,14 @@ def isentropic_youngs_modulus(self, direction):
         :rtype: float
         """
         direction = unit_normalize(direction)
-        Sijkl = self.full_isentropic_compliance_tensor
-        S = Sijkl.dot(direction).dot(direction).dot(direction).dot(direction)
+        S = np.einsum(
+            "ijkl, ...i, ...j, ...k, ...l",
+            self.full_isentropic_compliance_tensor,
+            direction,
+            direction,
+            direction,
+            direction,
+        )
         return 1.0 / S
 
     def isentropic_shear_modulus(self, plane_normal, shear_direction):
@@ -232,13 +239,16 @@ def isentropic_shear_modulus(self, plane_normal, shear_direction):
         assert (
             np.abs(plane_normal.dot(shear_direction)) < np.finfo(float).eps
         ), "plane_normal and shear_direction must be orthogonal"
-        Sijkl = self.full_isentropic_compliance_tensor
-        G = (
-            Sijkl.dot(shear_direction)
-            .dot(plane_normal)
-            .dot(shear_direction)
-            .dot(plane_normal)
+
+        G = np.einsum(
+            "ijkl, ...i, ...j, ...k, ...l",
+            self.full_isentropic_compliance_tensor,
+            shear_direction,
+            plane_normal,
+            shear_direction,
+            plane_normal,
         )
+
         return 0.25 / G
 
     def isentropic_poissons_ratio(self, axial_direction, lateral_direction):
@@ -250,14 +260,28 @@ def isentropic_poissons_ratio(self, axial_direction, lateral_direction):
 
         axial_direction = unit_normalize(axial_direction)
         lateral_direction = unit_normalize(lateral_direction)
-        assert (
-            np.abs(axial_direction.dot(lateral_direction)) < np.finfo(float).eps
-        ), "axial_direction and lateral_direction must be orthogonal"
 
-        Sijkl = self.full_isentropic_compliance_tensor
-        x = axial_direction
-        y = lateral_direction
-        nu = -(Sijkl.dot(y).dot(y).dot(x).dot(x) / Sijkl.dot(x).dot(x).dot(x).dot(x))
+        if len(axial_direction.shape) == 1:
+            assert (
+                np.abs(np.einsum("i, i", axial_direction, lateral_direction)) < 1.0e-12
+            ), "axial_direction and lateral_direction must be orthogonal"
+        else:
+            dot = np.einsum("...i, ...i->...", axial_direction, lateral_direction)
+            assert np.all(
+                np.abs(dot) < 1.0e-12
+            ), "axial_direction and lateral_direction must be orthogonal"
+
+        nu = -(
+            self.isentropic_youngs_modulus(axial_direction)
+            * np.einsum(
+                "ijkl, ...i, ...j, ...k, ...l",
+                self.full_isentropic_compliance_tensor,
+                lateral_direction,
+                lateral_direction,
+                axial_direction,
+                axial_direction,
+            )
+        )
         return nu
 
     def wave_velocities(self, propagation_direction):
@@ -267,15 +291,16 @@ def wave_velocities(self, propagation_direction):
         :rtype: list, containing the wave speeds and directions
             of particle motion relative to the stiffness tensor
         """
-        propagation_direction = unit_normalize(propagation_direction)
-
         Tik = self.christoffel_tensor(propagation_direction)
 
         eigenvalues, eigenvectors = np.linalg.eig(Tik)
 
-        idx = eigenvalues.argsort()[::-1]
-        eigenvalues = np.real(eigenvalues[idx])
-        eigenvectors = eigenvectors[:, idx]
+        # sort eigs by decreasing eigenvalue
+        idx = np.flip(eigenvalues.argsort(axis=-1), axis=-1)
+        vidx = np.expand_dims(idx, -1)
+        eigenvalues = np.take_along_axis(eigenvalues, idx, axis=-1)
+        eigenvectors = np.take_along_axis(eigenvectors, vidx, axis=-2)
+        eigenvalues = np.real(eigenvalues)
         velocities = np.sqrt(eigenvalues / self.density)
 
         return velocities, eigenvectors
diff --git a/tests/test_anisotropy.py b/tests/test_anisotropy.py
index 3bfc51bf..38e8860e 100644
--- a/tests/test_anisotropy.py
+++ b/tests/test_anisotropy.py
@@ -41,13 +41,14 @@ def test_isotropic(self):
         d1 = [1.0, 0.0, 0.0]
         d2 = [0.0, 1.0, 0.0]
 
-        beta_100 = m.isentropic_linear_compressibility(direction=d1)
-        E_100 = m.isentropic_youngs_modulus(direction=d1)
+        ds = np.array([d1, d2])
+        beta_100 = m.isentropic_linear_compressibility(direction=ds)[0]
+        E_100 = m.isentropic_youngs_modulus(direction=ds)[0]
         G_100_010 = m.isentropic_shear_modulus(plane_normal=d1, shear_direction=d2)
         nu_100_010 = m.isentropic_poissons_ratio(
             axial_direction=d1, lateral_direction=d2
         )
-        wave_speeds, wave_directions = m.wave_velocities(propagation_direction=d1)
+        wave_speeds, _ = m.wave_velocities(propagation_direction=d1)
         Vp, Vs1, Vs2 = wave_speeds
 
         array_2 = [