diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c
index febef46..8b837f9 100644
--- a/qpms/scatsystem.c
+++ b/qpms/scatsystem.c
@@ -22,6 +22,7 @@
 #include <pthread.h>
 #include "kahansum.h"
 #include "tolerances.h"
+#include "beyn.h"
 
 #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS
 #include "qpmsblas.h"
@@ -1897,3 +1898,33 @@ complex double *qpms_scatsys_scatter_solve(
   return f;
 }
 
+/// Wrapper for Beyn algorithm (non-periodic system)
+static int qpms_scatsys_finite_eval_Beyn_ImTS(complex double *target,
+    size_t m, complex double omega, void *params) {
+  const qpms_scatsys_t *ss = params;
+  qpms_scatsys_at_omega_t *ssw = qpms_scatsys_at_omega(ss, omega);
+  QPMS_ENSURE(ssw != NULL, "qpms_scatsys_at_omega() returned NULL");
+  QPMS_ENSURE(NULL != qpms_scatsysw_build_modeproblem_matrix_full(target, ssw),
+    "qpms_scatsysw_build_modeproblem_matrix_full() returned NULL");
+  qpms_scatsys_at_omega_free(ssw);
+  return QPMS_SUCCESS;
+}
+
+beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss,
+    complex double omega_centre, double omega_rr, double omega_ri, 
+    size_t contour_npoints, 
+    double rank_tol, size_t rank_sel_min, double res_tol) {
+
+  beyn_contour_t *contour = beyn_contour_ellipse(omega_centre,
+      omega_rr, omega_ri, contour_npoints);
+
+  beyn_result_t *result = beyn_solve(ss->fecv_size, 
+      ss->fecv_size /* possibly make smaller? */,
+      qpms_scatsys_finite_eval_Beyn_ImTS, NULL, (void *) ss,
+      contour, rank_tol, rank_sel_min, res_tol);
+  QPMS_ENSURE(result != NULL, "beyn_solve() returned NULL");
+
+  free(contour);
+  return result;
+}
+
diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h
index 49bc9b8..2aa1a99 100644
--- a/qpms/scatsystem.h
+++ b/qpms/scatsystem.h
@@ -343,31 +343,43 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
 		);
 
 /// Creates the full \f$ (I - TS) \f$ matrix of the scattering system.
+/**
+ * \returns \a target on success, NULL on error.
+ */
 complex double *qpms_scatsysw_build_modeproblem_matrix_full(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
 		complex double *target,
 		const qpms_scatsys_at_omega_t *ssw
 		);
 /// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form.
+/**
+ * \returns \a target on success, NULL on error.
+ */
 complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
 		complex double *target,
 		const qpms_scatsys_at_omega_t *ssw,
-		qpms_iri_t iri
+		qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
 		);
 /// Alternative implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed().
+/**
+ * \returns \a target on success, NULL on error.
+ */
 complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
 		complex double *target,
 		const qpms_scatsys_at_omega_t *ssw,
-		qpms_iri_t iri
+		qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
 		);
 /// Alternative (serial reference) implementation of qpms_scatsysw_build_modeproblem_matrix_irrep_packed().
+/**
+ * \returns \a target on success, NULL on error.
+ */
 complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
 		/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
 		complex double *target,
 		const qpms_scatsys_at_omega_t *ssw,
-		qpms_iri_t iri
+		qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
 		);
 
 /// LU factorisation (LAPACKE_zgetrf) result holder.
@@ -504,6 +516,31 @@ complex double *qpms_scatsysw_apply_Tmatrices_full(
 		const qpms_scatsys_at_omega_t *ssw
 		);
 
+struct beyn_result_t; // See beyn.h for full definition
+
+/// Searches for scattering system's eigenmodes using Beyn's algorithm.
+/**
+ * Currently, elliptical contour is used.
+ *
+ * TODO In the future, this will probably support irrep decomposition as well,
+ * but it does not make much sense for periodic / small systems, as in their
+ * case the bottleneck is the T-matrix and translation matrix evaluation
+ * rather than the linear algebra.
+ */
+
+struct beyn_result_t *qpms_scatsys_find_eigenmodes(
+		const qpms_scatsys_t *ss,
+		double eta, ///< Ewald sum parameter
+		const double *beta_, ///< k-vector of corresponding dimensionality, NULL/ignored for finite system.
+		complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
+		double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
+		double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
+		size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
+		double rank_tol, ///< (default: `1e-4`) TODO DOC.
+		size_t rank_sel_min, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
+		double res_tol ///< (default: `0.0`) TODO DOC.
+		);
+
 
 #if 0
 /// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points.