diff --git a/adcc/__init__.py b/adcc/__init__.py
index 9787c0d3..872963c0 100644
--- a/adcc/__init__.py
+++ b/adcc/__init__.py
@@ -75,7 +75,7 @@
 __url__ = "https://adc-connect.org"
 __authors__ = ["Michael F. Herbst", "Maximilian Scheurer", "Jonas Leitner",
                "Antonia Papapostolou", "Friederike Schneider",
-               "Adrian L. Dempwolff"]
+               "Adrian L. Dempwolff", "Adrian J. Müller"]
 __email__ = "developers@adc-connect.org"
 __contributors__ = []
 
diff --git a/adcc/solver/davidson.py b/adcc/solver/davidson.py
index d0588c7c..b3eca74d 100644
--- a/adcc/solver/davidson.py
+++ b/adcc/solver/davidson.py
@@ -78,7 +78,8 @@ def default_print(state, identifier, file=sys.stdout):
 def davidson_iterations(matrix, state, max_subspace, max_iter, n_ep, n_block,
                         is_converged, which, callback=None, preconditioner=None,
                         preconditioning_method="Davidson", debug_checks=False,
-                        residual_min_norm=None, explicit_symmetrisation=None):
+                        residual_min_norm=None, explicit_symmetrisation=None,
+                        max_subspace_iter=None):
     """Drive the davidson iterations
 
     Parameters
@@ -119,6 +120,8 @@ def davidson_iterations(matrix, state, max_subspace, max_iter, n_ep, n_block,
         Explicit symmetrisation to apply to new subspace vectors before
         adding them to the subspace. Allows to correct for loss of index
         or spin symmetries (type or instance)
+    max_subspace_iter : int, optional
+        Maximum number of iterations for diagonalizing the subspace matrix
     """
     if preconditioning_method not in ["Davidson", "Sleijpen-van-der-Vorst"]:
         raise ValueError("Only 'Davidson' and 'Sleijpen-van-der-Vorst' "
@@ -186,7 +189,8 @@ def callback(state, identifier):
                 # TODO Maybe play with precision a little here
                 # TODO Maybe use previous vectors somehow
                 v0 = None
-                rvals, rvecs = sla.eigsh(Ass, k=n_block, which=which, v0=v0)
+                rvals, rvecs = sla.eigsh(Ass, k=n_block, which=which, v0=v0,
+                                         maxiter=max_subspace_iter)
 
         with state.timer.record("residuals"):
             # Form residuals, A * SS * v - λ * SS * v = Ax * v + SS * (-λ*v)
@@ -350,7 +354,8 @@ def eigsh(matrix, guesses, n_ep=None, n_block=None, max_subspace=None,
           conv_tol=1e-9, which="SA", max_iter=70,
           callback=None, preconditioner=None,
           preconditioning_method="Davidson", debug_checks=False,
-          residual_min_norm=None, explicit_symmetrisation=IndexSymmetrisation):
+          residual_min_norm=None, explicit_symmetrisation=IndexSymmetrisation,
+          max_subspace_iter=None):
     """Davidson eigensolver for ADC problems
 
     Parameters
@@ -391,6 +396,12 @@ def eigsh(matrix, guesses, n_ep=None, n_block=None, max_subspace=None,
         Minimal norm a residual needs to have in order to be accepted as
         a new subspace vector
         (defaults to 2 * len(matrix) * machine_expsilon)
+    explicit_symmetrisation
+        Explicit symmetrisation to apply to new subspace vectors before
+        adding them to the subspace. Allows to correct for loss of index
+        or spin symmetries (type or instance)
+    max_subspace_iter : int, optional
+        Maximum number of iterations for diagonalizing the subspace matrix
     """
     if not isinstance(matrix, AdcMatrixlike):
         raise TypeError("matrix is not of type AdcMatrixlike")
@@ -452,7 +463,8 @@ def convergence_test(state):
                         preconditioning_method=preconditioning_method,
                         debug_checks=debug_checks,
                         residual_min_norm=residual_min_norm,
-                        explicit_symmetrisation=explicit_symmetrisation)
+                        explicit_symmetrisation=explicit_symmetrisation,
+                        max_subspace_iter=max_subspace_iter)
     return state
 
 
diff --git a/adcc/workflow.py b/adcc/workflow.py
index 6c79fb4c..a9d27bf6 100644
--- a/adcc/workflow.py
+++ b/adcc/workflow.py
@@ -139,6 +139,9 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None,
         Maximal subspace size
     max_iter : int, optional
         Maximal number of iterations
+    max_subspace_iter : int, optional
+        Maximum number of iterations for diagonalizing the subspace matrix
+        during the Davidson solver procedure
 
     Returns
     -------
diff --git a/pyproject.toml b/pyproject.toml
index cb748ef4..dcf2af0a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -15,6 +15,7 @@ authors = [
 	{name = "Antonia Papapostolou"},
 	{name = "Friederike Schneider"},
 	{name = "Adrian L. Dempwolff"},
+	{name = "Adrian J. Müller"},
 	{email = "developers@adc-connect.org"},
 ]
 license = "GPL-3.0-or-later"