diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx
index f7f853c..fed7007 100644
--- a/qpms/qpms_c.pyx
+++ b/qpms/qpms_c.pyx
@@ -589,6 +589,56 @@ cdef class ScatteringSystem:
                 self.s, qpms_incfield_planewave, <void *>&p, 0)
         return target_np
 
+    def find_modes(self, cdouble omega_centre, double omega_rr, double omega_ri,
+            size_t contour_points = 20, double rank_tol = 1e-4, size_t rank_min_sel=1,
+            double res_tol = 0):
+        """
+        Attempts to find the eigenvalues and eigenvectors using Beyn's algorithm.
+        """
+
+        cdef beyn_result_t *res = qpms_scatsys_finite_find_eigenmodes(self.s, 
+                omega_centre, omega_rr, omega_ri, contour_points, 
+                rank_tol, rank_min_sel, res_tol)
+        if res == NULL: raise RuntimeError
+
+        cdef size_t neig = res[0].neig
+        cdef size_t vlen = res[0].vlen # should be equal to self.s.fecv_size
+
+        cdef np.ndarray[complex, ndim=1] eigval = np.empty((neig,), dtype=complex)
+        cdef cdouble[::1] eigval_v = eigval
+        cdef np.ndarray[complex, ndim=1] eigval_err = np.empty((neig,), dtype=complex)
+        cdef cdouble[::1] eigval_err_v = eigval_err
+        cdef np.ndarray[double, ndim=1] residuals = np.empty((neig,), dtype=np.double)
+        cdef double[::1] residuals_v = residuals
+        cdef np.ndarray[complex, ndim=2] eigvec = np.empty((neig,vlen),dtype=complex)
+        cdef cdouble[:,::1] eigvec_v = eigvec
+        cdef np.ndarray[double, ndim=1] ranktest_SV = np.empty((vlen), dtype=np.double)
+        cdef double[::1] ranktest_SV_v = ranktest_SV
+
+        for i in range(neig):
+            eigval_v[i] = res[0].eigval[i]
+            eigval_err_v[i] = res[0].eigval_err[i]
+            residuals_v[i] = res[0].residuals[i]
+            for j in range(vlen):
+                eigvec_v[i,j] = res[0].eigvec[i*vlen + j]
+        for i in range(vlen):
+            ranktest_SV_v[i] = res[0].ranktest_SV[i]
+
+        zdist = eigval - omega_centre
+        eigval_inside_metric = np.hypot(zdist.real / omega_rr, zdist.imag / omega_ri)
+
+        beyn_result_free(res)
+        retdict = {
+                'eigval':eigval,
+                'eigval_inside_metric':eigval_inside_metric,
+                'eigvec':eigvec,
+                'residuals':residuals,
+                'eigval_err':eigval_err,
+                'ranktest_SV':ranktest_SV,
+        }
+
+        return retdict
+
 cdef class _ScatteringSystemAtOmega:
     '''
     Wrapper over the C qpms_scatsys_at_omega_t structure
diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd
index a0d892a..f408f5a 100644
--- a/qpms/qpms_cdefs.pxd
+++ b/qpms/qpms_cdefs.pxd
@@ -598,6 +598,9 @@ cdef extern from "scatsystem.h":
     cdouble *qpms_scatsys_scatter_solve(cdouble *target_f, const cdouble *a_inc, qpms_ss_LU ludata)
     const qpms_vswf_set_spec_t *qpms_ss_bspec_tmi(const qpms_scatsys_t *ss, qpms_ss_tmi_t tmi)
     const qpms_vswf_set_spec_t *qpms_ss_bspec_pi(const qpms_scatsys_t *ss, qpms_ss_pi_t pi)
+    beyn_result_t *qpms_scatsys_finite_find_eigenmodes(const qpms_scatsys_t *ss, cdouble omega_centre,
+            double omega_rr, double omega_ri, size_t contour_npoints, double rank_tol, size_t rank_sel_min,
+            double res_tol)
 
 cdef extern from "ewald.h":
     struct qpms_csf_result:
diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h
index 2aa1a99..54b9852 100644
--- a/qpms/scatsystem.h
+++ b/qpms/scatsystem.h
@@ -518,6 +518,27 @@ complex double *qpms_scatsysw_apply_Tmatrices_full(
 
 struct beyn_result_t; // See beyn.h for full definition
 
+/// Searches for finite 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_finite_find_eigenmodes(
+		const qpms_scatsys_t *ss,
+		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
 /// Searches for scattering system's eigenmodes using Beyn's algorithm.
 /**
  * Currently, elliptical contour is used.
@@ -527,7 +548,6 @@ struct beyn_result_t; // See beyn.h for full definition
  * 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
@@ -540,6 +560,7 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
 		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.
 		);
+#endif
 
 
 #if 0