diff --git a/qpms/beyn.c b/qpms/beyn.c
index fbe9763..3dce3a8 100644
--- a/qpms/beyn.c
+++ b/qpms/beyn.c
@@ -417,7 +417,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function,
    }
 
 */
-  return 0;
+  return K;
 }
 
 /***************************************************************/
diff --git a/qpms/beyn.h b/qpms/beyn.h
index 61ec3b4..6b9f69f 100644
--- a/qpms/beyn.h
+++ b/qpms/beyn.h
@@ -33,6 +33,7 @@
 
 #include <complex.h>
 #include <gsl/gsl_matrix.h>
+#include <gsl/gsl_complex_math.h>
 
 /// User-supplied function that provides the operator M(z) whose "roots" are to be found.
 typedef int (*beyn_function_M_t)(gsl_matrix_complex *target_M, complex double z, void *params);
@@ -83,4 +84,8 @@ int BeynSolve(BeynSolver *Solver,
               beyn_function_M_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *params,
               double complex z0, double Rx, double Ry, int N);
 
+
+static inline complex double gsl_complex_tostd(gsl_complex z) { return GSL_REAL(z) + I*GSL_IMAG(z); }
+static inline gsl_complex gsl_complex_fromstd(complex double z) { return gsl_complex_rect(creal(z), cimag(z)); }
+
 #endif // BEYN_H
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 6e6797b..eb0b17e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -15,7 +15,11 @@ add_executable( test_scalar_ewald32 test_scalar_ewald32.c )
 target_link_libraries( test_scalar_ewald32 qpms gsl lapacke amos m )
 target_include_directories( test_scalar_ewald32 PRIVATE .. )
 
-add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array )
+add_executable( tbeyn tbeyn.c )
+target_link_libraries( tbeyn qpms gsl lapacke amos m )
+target_include_directories( tbeyn PRIVATE .. )
+
+add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array tbeyn )
 
 add_test( NAME single_vs_array_translation_coeffs COMMAND test_single_translations_vs_calc )
 
diff --git a/tests/tbeyn.c b/tests/tbeyn.c
new file mode 100644
index 0000000..a78113f
--- /dev/null
+++ b/tests/tbeyn.c
@@ -0,0 +1,39 @@
+#include <qpms/beyn.h>
+#include <stdio.h>
+
+
+// Matrix as in Beyn, section 4.11
+int M_function(gsl_matrix_complex *target, complex double z, void *no_params) {
+  int m = target->size1;
+
+  gsl_complex d = gsl_complex_fromstd(  2*m - 4*z / (6*m) );
+  gsl_complex od =  gsl_complex_fromstd( -m - z / (6*m) );
+
+  gsl_matrix_complex_set_zero(target);
+  for (int i = 0; i < m; ++i) {
+    gsl_matrix_complex_set(target, i, i, (gsl_complex) d);
+    if(i > 0) gsl_matrix_complex_set(target, i, i-1, (gsl_complex) od);
+    if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, (gsl_complex) od);
+  }
+
+  return 0;
+}
+
+int main() {
+  complex double z0 = 150+2*I;
+  double Rx = 148, Ry = 148;
+  int L = 10, N = 50, dim = 400;
+  BeynSolver * solver = CreateBeynSolver(dim, L);
+
+  int K = BeynSolve(solver, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
+      z0, Rx, Ry, N);
+  printf("Found %d eigenvalues:\n", K);
+  for (int i = 0; i < K; ++i) {
+    gsl_complex eig = gsl_vector_complex_get(solver->Eigenvalues, i);
+    printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
+  }
+  DestroyBeynSolver(solver);
+  return 0;
+}
+
+