diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt
index 575aca8..2a63a3f 100644
--- a/qpms/CMakeLists.txt
+++ b/qpms/CMakeLists.txt
@@ -12,7 +12,7 @@ add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c e
 	ewaldsf.c pointgroups.c latticegens.c
 	lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
 	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
-	lll.c
+	lll.c beyn.c
 	)
 use_c99()
 
diff --git a/qpms/beyn.c b/qpms/beyn.c
new file mode 100644
index 0000000..28911be
--- /dev/null
+++ b/qpms/beyn.c
@@ -0,0 +1,667 @@
+#define STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1]
+
+#define cg2s(x) gsl_complex_tostd(x)
+#define cs2g(x) gsl_complex_fromstd(x)
+
+#include <complex.h>
+#include <lapacke.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "qpms_error.h"
+
+// Maybe GSL works?
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_eigen.h>
+
+#include "beyn.h"
+#define SQ(x) ((x)*(x))
+
+STATIC_ASSERT((sizeof(lapack_complex_double) == sizeof(gsl_complex)), lapack_and_gsl_complex_must_be_consistent);
+
+
+typedef struct BeynSolver
+{
+  int M;   // dimension of matrices
+  int L;   // number of columns of VHat matrix
+
+  gsl_vector_complex *eigenvalues, *eigenvalue_errors;
+  gsl_matrix_complex *eigenvectors;
+  gsl_matrix_complex *A0, *A1, *A0_coarse, *A1_coarse, *MInvVHat;
+  gsl_matrix_complex *VHat;
+  gsl_vector *Sigma, *residuals;
+} BeynSolver;
+
+// constructor, destructor
+BeynSolver *BeynSolver_create(int M, int L);
+void BeynSolver_free(BeynSolver *solver);
+
+// reset the random matrix VHat used in Beyn's algorithm
+void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed);
+
+// Uniformly random number from interval [a, b].
+static double randU(double a, double b) { return a + (b-a) * random() * (1. / RAND_MAX); }
+
+// Random number from normal distribution (via Box-Muller transform, which is enough for our purposes).
+static double randN(double Sigma, double Mu) {
+  double u1 = randU(0, 1);
+  double u2 = randU(0, 1);
+  return Mu + Sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2);
+}
+
+static complex double zrandN(double sigma, double mu){
+  return randN(sigma, mu) + I*randN(sigma, mu);
+}
+
+static inline double dsq(double a) { return a * a; }
+
+static _Bool beyn_contour_ellipse_inside_test(struct beyn_contour_t *c, complex double z) {
+  double rRe = c->z_dz[c->n][0];
+  double rIm = c->z_dz[c->n][1];
+  complex double zn = z - c->centre;
+  return dsq(creal(zn)/rRe) + dsq(cimag(zn)/rIm) < 1;
+}
+
+beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double rIm, size_t n) 
+{
+  beyn_contour_t *c;
+  QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0]));
+  c->centre = centre;
+  c->n = n;
+  for(size_t i = 0; i < n; ++i) {
+    double t = i*2*M_PI/n;
+    double st = sin(t), ct = cos(t);
+    c->z_dz[i][0] = centre + ct*rRe + I*st*rIm;
+    c->z_dz[i][1] = (-rRe*st + I*rIm*ct) * (2*M_PI / n);
+  }
+  // We hide the half-axes metadata after the discretisation points.
+  c->z_dz[n][0] = rRe;
+  c->z_dz[n][1] = rIm;
+  c->inside_test = beyn_contour_ellipse_inside_test;
+  return c;
+}
+
+
+// Sets correct sign to zero for a given branch cut orientation
+static inline complex double 
+align_zero(complex double z, beyn_contour_halfellipse_orientation or) 
+{
+  // Maybe redundant, TODO check the standard.
+  const double positive_zero = copysign(0., +1.);
+  const double negative_zero = copysign(0., -1.);
+  switch(or) { // ensure correct zero signs; CHECKME!!!
+    case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS:
+      if(creal(z) == 0 && signbit(creal(z)))
+        z = positive_zero + I * cimag(z);
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS:
+      if(creal(z) == 0 && !signbit(creal(z)))
+        z = negative_zero + I * cimag(z);
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS:
+      if(cimag(z) == 0 && signbit(cimag(z)))
+        z = creal(z) + I * positive_zero;
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS:
+      if(cimag(z) == 0 && !signbit(cimag(z)))
+        z = creal(z) + I * negative_zero;
+      break;
+    default: QPMS_WTF;
+  }
+  return z;
+}
+
+
+beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe,
+    double rIm, size_t n, beyn_contour_halfellipse_orientation or)
+{
+  beyn_contour_t *c;
+  QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0])
+      + sizeof(beyn_contour_halfellipse_orientation));
+  c->centre = centre;
+  c->n = n;
+  const size_t nline = n/2;
+  const size_t narc = n - nline;
+  complex double faktor;
+  double l = rRe, h = rIm;
+  switch(or) {
+    case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS:
+      faktor = -I;
+      l = rIm; h = rRe;
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS:
+      faktor = I;
+      l = rIm; h = rRe;
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS:
+      faktor = 1;
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS:
+      faktor = -1;
+      break;
+    default: QPMS_WTF;
+  }
+
+  for(size_t i = 0; i < narc; ++i) {
+    double t = (i+0.5)*M_PI/narc;
+    double st = sin(t), ct = cos(t);
+    c->z_dz[i][0] = centre + faktor*(ct*l + I*st*h);
+    c->z_dz[i][1] = faktor * (-l*st + I*h*ct) * (M_PI / narc);
+  }
+  for(size_t i = 0; i < nline; ++i) {
+    double t = 0.5 * (1 - (double) nline) + i;
+    c->z_dz[narc + i][0] = align_zero(centre + faktor * t * 2 * l / nline, or);
+    c->z_dz[narc + i][1] = faktor * 2 * l / nline;
+  }
+  // We hide the half-axes metadata after the discretisation points.
+  c->z_dz[n][0] = rRe;
+  c->z_dz[n][1] = rIm;
+  // ugly...
+  *((beyn_contour_halfellipse_orientation *) &c->z_dz[n+1][0]) = or;
+  c->inside_test = NULL; // TODO beyn_contour_halfellipse_inside_test;
+  return c;
+}
+
+beyn_contour_t *beyn_contour_kidney(complex double centre, double rRe,
+    double rIm, const double rounding, const size_t n, beyn_contour_halfellipse_orientation or)
+{
+  beyn_contour_t *c;
+  QPMS_ENSURE(rounding >= 0 && rounding < .5, "rounding must lie in the interval [0, 0.5)");
+  QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0])
+      + sizeof(beyn_contour_halfellipse_orientation));
+  c->centre = centre;
+  c->n = n;
+  complex double faktor;
+  double l = rRe, h = rIm;
+  switch(or) {
+    case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS:
+      faktor = -I;
+      l = rIm; h = rRe;
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_RE_MINUS:
+      faktor = I;
+      l = rIm; h = rRe;
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_IM_PLUS:
+      faktor = 1;
+      break;
+    case BEYN_CONTOUR_HALFELLIPSE_IM_MINUS:
+      faktor = -1;
+      break;
+    default: QPMS_WTF;
+  }
+
+  // Small circle centre coordinates.
+  const double y = rounding * h; // distance from the cut / straight line
+  const double x = sqrt(SQ(h - y) - SQ(y)); 
+
+  const double alpha = asin(y/(h-y));
+  const double ar = l/h; // aspect ratio
+
+  // Parameter range (equal to the contour length if ar == 1)
+  const double tmax = 2 * (x + (M_PI_2 + alpha) * y + (M_PI_2 - alpha) * h);
+  const double dt = tmax / n;
+
+  size_t i = 0;
+  double t;
+  // Straight line, first part
+  double t_lo = 0, t_hi = x;
+  for(; t = i * dt, t <= t_hi; ++i) { 
+    c->z_dz[i][0] = align_zero(centre + (t - t_lo) * ar * faktor, or);
+    c->z_dz[i][1] = dt * ar * faktor;
+  }
+  // First small arc
+  t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y;
+  for(; t = i * dt, t < t_hi; ++i) {
+    double phi = (t - t_lo) / y - M_PI_2;
+    c->z_dz[i][0] = centre + (ar * (x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor;
+    c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
+  }
+  // Big arc
+  t_lo = t_hi; t_hi = t_lo + (M_PI - 2 * alpha) * h;
+  for(; t = i * dt, t < t_hi; ++i) {
+    double phi = (t - t_lo) / h + alpha;
+    c->z_dz[i][0] = centre + (ar * (h * cos(phi)) + h * sin(phi) * I) * faktor;
+    c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
+  }
+  // Second small arc
+  t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y;
+  for(; t = i * dt, t < t_hi; ++i) {
+    double phi = (t - t_lo) / y + M_PI - alpha;
+    c->z_dz[i][0] = centre + (ar * (- x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor;
+    c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
+  }
+  // Straight line, second part
+  t_lo = t_hi; t_hi = tmax;
+  for(; t = i * dt, i < n; ++i) {
+    c->z_dz[i][0] = align_zero(centre + (t - t_lo - x) * ar * faktor, or);
+    c->z_dz[i][1] = dt * ar * faktor;
+  }
+  
+#if 0 // TODO later
+  // We hide the half-axes metadata after the discretisation points.
+  c->z_dz[n][0] = rRe;
+  c->z_dz[n][1] = rIm;
+  // ugly...
+  *((beyn_contour_halfellipse_orientation *) &c->z_dz[n+1][0]) = or;
+#endif
+  c->inside_test = NULL; // TODO beyn_contour_halfellipse_inside_test;
+  return c;
+}
+
+void beyn_result_gsl_free(beyn_result_gsl_t *r) {
+  if(r) {
+    gsl_vector_complex_free(r->eigval);
+    gsl_vector_complex_free(r->eigval_err);
+    gsl_vector_free(r->residuals);
+    gsl_matrix_complex_free(r->eigvec);
+    gsl_vector_free(r->ranktest_SV);
+    free(r);
+  }
+}
+
+BeynSolver *BeynSolver_create(int M, int L)
+{
+  BeynSolver *solver= (BeynSolver *)malloc(sizeof(*solver));
+
+  solver->M = M;
+  solver->L = L;
+  QPMS_ENSURE(L <= M, "We expect L <= M, but we got L = %d, M = %d ", L, M);
+
+  // storage for eigenvalues and eigenvectors
+  solver->eigenvalues  = gsl_vector_complex_calloc(L);
+  solver->eigenvalue_errors     = gsl_vector_complex_calloc(L);
+  solver->residuals    = gsl_vector_calloc(L);
+  solver->eigenvectors = gsl_matrix_complex_calloc(M, L);
+
+  // storage for singular values, random VHat matrix, etc. used in algorithm
+  solver->A0           = gsl_matrix_complex_calloc(M,L);
+  solver->A1           = gsl_matrix_complex_calloc(M,L);
+  solver->A0_coarse     = gsl_matrix_complex_calloc(M,L);
+  solver->A1_coarse     = gsl_matrix_complex_calloc(M,L);
+  solver->MInvVHat     = gsl_matrix_complex_calloc(M,L);
+  solver->VHat         = gsl_matrix_complex_calloc(M,L);
+  solver->Sigma        = gsl_vector_calloc(L);
+  // Beyn Step 1: Generate random matrix VHat
+  BeynSolver_srandom(solver,(unsigned)time(NULL));
+
+  return solver;
+
+}
+
+void BeynSolver_free(BeynSolver *solver)
+{
+  gsl_vector_complex_free(solver->eigenvalues);
+  gsl_vector_complex_free(solver->eigenvalue_errors);  
+  gsl_matrix_complex_free(solver->eigenvectors);
+
+  gsl_matrix_complex_free(solver->A0);
+  gsl_matrix_complex_free(solver->A1);
+  gsl_matrix_complex_free(solver->A0_coarse);
+  gsl_matrix_complex_free(solver->A1_coarse);
+  gsl_matrix_complex_free(solver->MInvVHat);
+  gsl_vector_free(solver->Sigma);
+  gsl_vector_free(solver->residuals);
+  gsl_matrix_complex_free(solver->VHat);
+
+  free(solver);
+}
+
+void BeynSolver_free_partial(BeynSolver *solver) 
+{
+  gsl_matrix_complex_free(solver->A0);
+  gsl_matrix_complex_free(solver->A1);
+  gsl_matrix_complex_free(solver->A0_coarse);
+  gsl_matrix_complex_free(solver->A1_coarse);
+  gsl_matrix_complex_free(solver->MInvVHat);
+  gsl_matrix_complex_free(solver->VHat);
+
+  free(solver);
+}
+
+void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed)
+{
+  if (RandSeed==0) 
+    RandSeed=time(0);
+  srandom(RandSeed);
+  gsl_matrix_complex *VHat=solver->VHat;
+  for(int nr=0; nr<VHat->size1; nr++)
+    for(int nc=0; nc<VHat->size2; nc++)
+      gsl_matrix_complex_set(VHat,nr,nc,cs2g(zrandN(1, 0)));
+
+}
+
+
+/*
+ * linear-algebra manipulations on the A0 and A1 matrices 
+ * (obtained via numerical quadrature) to extract eigenvalues 
+ * and eigenvectors                                
+ */
+
+static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_function, 
+    void *Params,
+    gsl_matrix_complex *A0, gsl_matrix_complex *A1, double complex z0,
+    gsl_vector_complex *eigenvalues, gsl_matrix_complex *eigenvectors, const double rank_tol, const double res_tol)
+{
+  const size_t m = solver->M;
+  const size_t l = solver->L;
+  gsl_vector *Sigma = solver->Sigma;
+
+  int verbose = 1; // TODO
+
+  // Beyn Step 3: Compute SVD: A0 = V0_full * Sigma * W0T_full
+  if(verbose) printf(" Beyn: computing SVD...\n");
+  gsl_matrix_complex* V0_full = gsl_matrix_complex_alloc(m,l);
+  gsl_matrix_complex_memcpy(V0_full,A0);
+  gsl_matrix_complex* W0T_full = gsl_matrix_complex_alloc(l,l);
+  QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
+  QPMS_ENSURE(V0_full->size1 >= V0_full->size2, "m = %zd, l = %zd, what the hell?");
+  QPMS_ENSURE_SUCCESS(LAPACKE_zgesdd(LAPACK_ROW_MAJOR, // A = U*Σ*conjg(V') 
+        'O' /*jobz, 'O' overwrites a with U and saves conjg(V') in vt if m >= n, i.e. if M >= L, which holds*/,
+        V0_full->size1 /* m, number of rows */,
+        V0_full->size2 /* n, number of columns */,
+        (lapack_complex_double *)(V0_full->data) /*a*/,
+        V0_full->tda /*lda*/,
+        Sigma->data /*s*/, 
+        NULL /*u; not used*/, 
+        m /*ldu; should not be really used but lapacke requires it for some obscure reason*/,
+        (lapack_complex_double *)W0T_full->data /*vt*/,
+        W0T_full->tda /*ldvt*/
+        ));
+
+
+  // Beyn Step 4: Rank test for Sigma
+  // compute effective rank K (number of eigenvalue candidates)
+  int K=0;
+  for (int k=0; k<Sigma->size /* this is l, actually */; k++) { 
+    if (verbose) printf("Beyn: SV(%d)=%e\n",k,gsl_vector_get(Sigma, k));
+    if (gsl_vector_get(Sigma, k) > rank_tol)
+      K++;
+  }
+  if (verbose)printf(" Beyn: %d/%zd relevant singular values\n",K,l);
+  if (K==0) { 
+    QPMS_WARN("no singular values found in Beyn eigensolver\n");
+    return 0;
+  }
+
+  // Beyn step 5:  B = V0' * A1 * W0 * Sigma^-1
+  // set V0, W0T = matrices of first K right, left singular vectors
+  gsl_matrix_complex *V0 = gsl_matrix_complex_alloc(m,K);
+  gsl_matrix_complex *W0T= gsl_matrix_complex_alloc(K,l);
+
+  for (int k = 0; k < K; ++k) {
+    gsl_vector_complex_view tmp;
+    tmp = gsl_matrix_complex_column(V0_full, k);
+    gsl_matrix_complex_set_col(V0, k, &(tmp.vector));
+    tmp = gsl_matrix_complex_row(W0T_full, k);
+    gsl_matrix_complex_set_row(W0T, k, &(tmp.vector));
+  }
+
+  gsl_matrix_complex_free(V0_full);
+  gsl_matrix_complex_free(W0T_full);
+
+  gsl_matrix_complex *TM2 = gsl_matrix_complex_calloc(K,l);
+  gsl_matrix_complex *B = gsl_matrix_complex_calloc(K,K);
+
+  if(verbose) printf(" Multiplying V0*A1->TM...\n");
+
+  const gsl_complex one = gsl_complex_rect(1,0);
+  const gsl_complex zero = gsl_complex_rect(0,0);
+  gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one,
+      V0, A1, zero, TM2);
+
+  if(verbose) printf(" Multiplying TM*W0T->B...\n");
+
+  gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one,
+      TM2, W0T, zero, B);
+
+  gsl_matrix_complex_free(W0T);
+  gsl_matrix_complex_free(TM2);
+
+  if(verbose) printf(" Scaling B <- B*Sigma^{-1}\n");
+  gsl_vector_complex *tmp = gsl_vector_complex_calloc(K);
+  for(int i = 0; i < K; i++) {
+    gsl_matrix_complex_get_col(tmp, B, i);
+    gsl_vector_complex_scale(tmp, gsl_complex_rect(1.0/gsl_vector_get(Sigma,i), 0));
+    gsl_matrix_complex_set_col(B,i,tmp);
+  }
+  gsl_vector_complex_free(tmp);
+
+  //for(int m=0; m<K; m++)                  //  B <- B * Sigma^{-1}
+  
+  // Beyn step 6.
+  // Eigenvalue decomposition B -> S*Lambda*S'
+  /* According to Beyn's paper (Algorithm 1), one should check conditioning
+   * of the eigenvalues; if they are ill-conditioned, one should perform
+   * a procedure involving Schur decomposition and reordering.
+   *
+   * Beyn refers to MATLAB routines eig, condeig, schur, ordschur.
+   * I am not sure about the equivalents in LAPACK, TODO check zgeevx, zgeesx.
+   */
+  if(verbose) printf(" Eigensolving (%i,%i)\n",K,K);
+
+  gsl_vector_complex *Lambda = gsl_vector_complex_alloc(K); // eigenvalues
+  gsl_matrix_complex *S = gsl_matrix_complex_alloc(K,K); // eigenvectors
+
+  QPMS_ENSURE(Sigma->stride == 1, "Sigma vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
+  QPMS_ENSURE(Lambda->stride == 1, "Lambda vector stride must be 1 for LAPACKE_zgesdd, got %zd.", Sigma->stride);
+  QPMS_ENSURE_SUCCESS(LAPACKE_zgeev(
+        LAPACK_ROW_MAJOR,
+        'N' /* jobvl; don't compute left eigenvectors */,
+        'V' /* jobvr; do compute right eigenvectors */,
+        K   /* n */,
+        (lapack_complex_double *)(B->data) /* a */,
+        B->tda  /* lda */,
+        (lapack_complex_double *) Lambda->data /* w */,
+        NULL /* vl */,
+        m /* ldvl, not used by for some reason needed */,
+        (lapack_complex_double *)(S->data)/* vr */,
+        S->tda/* ldvr */
+        ));
+
+  gsl_matrix_complex_free(B);
+
+  // V0S <- V0*S
+  printf("Multiplying V0*S...\n");
+  gsl_matrix_complex *V0S = gsl_matrix_complex_alloc(m, K);
+  QPMS_ENSURE_SUCCESS(gsl_blas_zgemm(CblasNoTrans, CblasNoTrans,
+        one, V0, S, zero, V0S));
+
+  gsl_matrix_complex_free(V0);
+
+  // FIXME!!! make clear relation between KRetained and K in the results!
+  // (If they differ, there are possibly some spurious eigenvalues.)
+  int KRetained = 0;
+  gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m);
+  gsl_vector_complex *MVk = gsl_vector_complex_alloc(m);
+  for (int k = 0; k < K; ++k) {
+    const gsl_complex zgsl = gsl_complex_add(gsl_complex_rect(creal(z0), cimag(z0)), gsl_vector_complex_get(Lambda, k));
+    const complex double z = GSL_REAL(zgsl) + I*GSL_IMAG(zgsl);
+    gsl_vector_complex_const_view Vk = gsl_matrix_complex_const_column(V0S, k);
+
+    double residual = 0;
+    if(res_tol > 0) {
+      QPMS_ENSURE_SUCCESS(M_function(Mmat, z, Params));
+      QPMS_ENSURE_SUCCESS(gsl_blas_zgemv(CblasNoTrans, one, Mmat, &(Vk.vector), zero, MVk));
+      residual = gsl_blas_dznrm2(MVk);
+      if (verbose) printf("Beyn: Residual(%i)=%e\n",k,residual);
+    }
+    if (res_tol > 0 && residual > res_tol) continue;
+
+    gsl_vector_complex_set(eigenvalues, KRetained, zgsl);
+    if(eigenvectors) {
+      gsl_matrix_complex_set_col(eigenvectors, KRetained, &(Vk.vector));
+      gsl_vector_set(solver->residuals, KRetained, residual);
+    }
+    ++KRetained;
+  }
+
+  gsl_matrix_complex_free(V0S);
+  gsl_matrix_complex_free(Mmat);
+  gsl_vector_complex_free(MVk);
+  gsl_matrix_complex_free(S);
+  gsl_vector_complex_free(Lambda);
+
+  return KRetained;
+}
+
+
+beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l,
+    beyn_function_M_gsl_t M_function, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat_function,
+    void *params, const beyn_contour_t *contour,
+    double rank_tol, double res_tol)
+{  
+  BeynSolver *solver = BeynSolver_create(m, l);
+
+  gsl_matrix_complex *A0           = solver->A0;
+  gsl_matrix_complex *A1           = solver->A1;
+  gsl_matrix_complex *A0_coarse     = solver->A0_coarse;
+  gsl_matrix_complex *A1_coarse     = solver->A1_coarse;
+  gsl_matrix_complex *MInvVHat     = solver->MInvVHat;
+  gsl_matrix_complex *VHat         = solver->VHat;
+
+  /***************************************************************/
+  /* evaluate contour integrals by numerical quadrature to get   */
+  /* A0 and A1 matrices                                          */
+  /***************************************************************/
+
+  gsl_matrix_complex_set_zero(A0);
+  gsl_matrix_complex_set_zero(A1);
+  gsl_matrix_complex_set_zero(A0_coarse);
+  gsl_matrix_complex_set_zero(A1_coarse);
+  const size_t N = contour->n;
+  if(N & 1) QPMS_WARN("Contour discretisation point number is odd (%zd),"
+     " the error estimates might be a bit off.", N);
+
+
+  // Beyn Step 2: Computa A0, A1 
+  const complex double z0 = contour->centre;
+  for(int n=0; n<N; n++) { 
+    const complex double z = contour->z_dz[n][0];
+    const complex double dz = contour->z_dz[n][1];
+
+    gsl_matrix_complex_memcpy(MInvVHat, VHat);
+
+    if(M_inv_Vhat_function) {
+      QPMS_ENSURE_SUCCESS(M_inv_Vhat_function(MInvVHat, VHat, z, params));
+    } else {
+      lapack_int *pivot;
+      gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m,m);
+      QPMS_ENSURE_SUCCESS(M_function(Mmat, z, params));
+      QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * m);
+      QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR,
+            m /*m*/, m /*n*/,(lapack_complex_double *) Mmat->data /*a*/, Mmat->tda /*lda*/, pivot /*ipiv*/));
+      QPMS_ENSURE(MInvVHat->tda == l, "wut?");
+      QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N' /*trans*/,
+            m /*n*/, l/*nrhs*/, (lapack_complex_double *)(Mmat->data) /*a*/, Mmat->tda /*lda*/, pivot/*ipiv*/, 
+            (lapack_complex_double *)(MInvVHat->data) /*b*/, MInvVHat->tda/*ldb*/));
+
+      free(pivot);
+      gsl_matrix_complex_free(Mmat);
+    }
+
+    gsl_matrix_complex_scale(MInvVHat, cs2g(dz));
+    gsl_matrix_complex_add(A0, MInvVHat);
+    if((n%2)==0) {
+      gsl_matrix_complex_add(A0_coarse, MInvVHat);
+      gsl_matrix_complex_add(A0_coarse, MInvVHat);
+    }
+
+    gsl_matrix_complex_scale(MInvVHat, cs2g(z - z0)); // A_1 scaling as in Beyn's Remark 3.2(b) for numerical stability.
+    gsl_matrix_complex_add(A1, MInvVHat);
+    if((n%2)==0) {
+      gsl_matrix_complex_add(A1_coarse, MInvVHat);
+      gsl_matrix_complex_add(A1_coarse, MInvVHat);
+    }
+  }
+
+  gsl_vector_complex *eigenvalues  = solver->eigenvalues;
+  gsl_vector_complex *eigenvalue_errors     = solver->eigenvalue_errors;
+  gsl_matrix_complex *eigenvectors = solver->eigenvectors;
+
+  // Beyn Steps 3–6
+  int K = beyn_process_matrices(solver, M_function, params, A0, A1, z0, eigenvalues, eigenvectors, rank_tol, res_tol);
+  // Repeat Steps 3–6 with rougher contour approximation to get an error estimate.
+  int K_coarse = beyn_process_matrices(solver, M_function, params, A0_coarse, A1_coarse, z0, eigenvalue_errors, eigenvectors, rank_tol, res_tol);
+
+  gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors);
+  // Reid did also fabs on the complex and real parts ^^^.
+
+  beyn_result_gsl_t *result;
+  QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t));
+  result->eigval = solver->eigenvalues;
+  result->eigval_err = solver->eigenvalue_errors;
+  result->residuals = solver->residuals;
+  result->eigvec = solver->eigenvectors;
+  result->ranktest_SV = solver->Sigma;
+  result->neig = K;
+
+  BeynSolver_free_partial(solver);
+
+  return result;
+}
+
+// Wrapper of pure C array M-matrix function to GSL.
+
+struct beyn_function_M_carr2gsl_param {
+  beyn_function_M_t M_function;
+  beyn_function_M_inv_Vhat_t M_inv_Vhat_function;
+  void *param;
+};
+
+static int beyn_function_M_carr2gsl(gsl_matrix_complex *target_M, complex double z, void *params) 
+{
+  struct beyn_function_M_carr2gsl_param *p = params;
+  // These could rather be asserts.
+  QPMS_ENSURE(target_M->size2 == target_M->tda, "Target GSL matrix is not a C-contiguous array. This is a bug, please report!");
+  QPMS_ENSURE(target_M->size1 == target_M->size2, "Target is not a square matrix. This is a bug, please report!");
+  return p->M_function((complex double *) target_M->data, target_M->size1, z, p->param);
+}
+
+static int beyn_function_M_inv_Vhat_carr2gsl(gsl_matrix_complex *target,
+    const gsl_matrix_complex *Vhat, complex double z, void *params) 
+{
+  QPMS_UNTESTED;
+  struct beyn_function_M_carr2gsl_param *p = params;
+  // These could rather be asserts.
+  QPMS_ENSURE(target->size2 == target->tda, "Target GSL matrix is not a C-contiguous array. This is a bug, please report!");
+  QPMS_ENSURE(Vhat->size2 == Vhat->tda, "Source GSL matrix is not a C-contiguous array. This is a bug, please report!");
+  // TODO here I could also check whether Vhat and target have compatible sizes.
+  return p->M_inv_Vhat_function((complex double *) target->data, target->size1, target->size2, 
+      (complex double *) Vhat->data, z, p->param);
+}
+
+beyn_result_t *beyn_solve(size_t m, size_t l, beyn_function_M_t M, beyn_function_M_inv_Vhat_t M_inv_Vhat,
+    void *params, const beyn_contour_t *contour, double rank_tol, double res_tol) {
+  struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params};
+  return beyn_result_from_beyn_result_gsl(
+    beyn_solve_gsl(m, l, beyn_function_M_carr2gsl,
+      (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, 
+      (void *) &p, contour, rank_tol, res_tol)
+    );
+}
+
+beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) {
+  struct beyn_result_t *result;
+  QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_t));
+  result->gsl = g;
+  result->neig = result->gsl->neig;
+  result->vlen = result->gsl->eigvec->size2;
+  result->eigval = (complex double *) result->gsl->eigval->data;
+  result->eigval_err = (complex double *) result->gsl->eigval_err->data;
+  result->residuals = result->gsl->residuals->data;
+  result->eigvec = (complex double *) result->gsl->eigvec->data;
+  result->ranktest_SV = result->gsl->ranktest_SV->data;
+  return result;
+}
+  
+void beyn_result_free(beyn_result_t *result) {
+  if(result)
+    beyn_result_gsl_free(result->gsl);
+  free(result);
+}
+
diff --git a/qpms/beyn.h b/qpms/beyn.h
new file mode 100644
index 0000000..25bc1c7
--- /dev/null
+++ b/qpms/beyn.h
@@ -0,0 +1,129 @@
+/** \file beyn.h
+ * \brief Beyn's algorithm for nonlinear eigenvalue problems.
+ */
+#ifndef BEYN_H
+#define BEYN_H
+
+#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.
+/** GSL matrix version */
+typedef int (*beyn_function_M_gsl_t)(gsl_matrix_complex *target_M, complex double z, void *params);
+
+/// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
+/** GSL matrix version */
+typedef int (*beyn_function_M_inv_Vhat_gsl_t)(gsl_matrix_complex *target_M_inv_Vhat,
+	       	const gsl_matrix_complex *Vhat, complex double z, void *params); 
+
+/// User-supplied function that provides the operator M(z) whose "roots" are to be found.
+/** Pure C array version */
+typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params);
+
+/// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
+/** Pure C array version */
+typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l,
+	       	const complex double *Vhat, complex double z, void *params); 
+
+/// Complex plane integration contour structure.
+typedef struct beyn_contour_t {
+	size_t n; ///< Number of discretisation points.
+	/// "Centre" of the contour.
+	/** 
+	 * This point is used in the rescaling of the \f$ A_1 \f$ matrix as in
+	 * Beyn's Remark 3.2 (b) in order to improve the numerical stability.
+	 * It does not have to be a centre in some strictly defined sense,
+	 * but it should be "somewhere around" where the contour is.
+	 */
+	complex double centre;
+	/// Function testing that a point \a z lies inside the contour (optional).
+	_Bool (*inside_test)(struct beyn_contour_t *, complex double z);
+	complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
+} beyn_contour_t;
+
+/// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes.
+/** Free using free(). */
+beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints);
+
+
+typedef enum {
+	BEYN_CONTOUR_HALFELLIPSE_RE_PLUS = 3,
+	BEYN_CONTOUR_HALFELLIPSE_RE_MINUS = 1,
+	BEYN_CONTOUR_HALFELLIPSE_IM_PLUS = 0,
+	BEYN_CONTOUR_HALFELLIPSE_IM_MINUS = 2,
+} beyn_contour_halfellipse_orientation;
+
+
+/// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes.
+/** Free using free(). */
+beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints,
+		beyn_contour_halfellipse_orientation or);
+
+/// Similar to halfellipse but with rounded corners.
+beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im, 
+	double rounding, ///< Must be in interval [0, 0.5)
+	size_t n, beyn_contour_halfellipse_orientation or);
+
+
+/// Beyn algorithm result structure (GSL matrix/vector version).
+typedef struct beyn_result_gsl_t {
+	size_t neig; ///< Number of eigenvalues found (a bit redundant?).
+	gsl_vector_complex *eigval;
+	gsl_vector_complex *eigval_err;
+	gsl_vector *residuals;
+	gsl_matrix_complex *eigvec;
+	gsl_vector *ranktest_SV;
+} beyn_result_gsl_t;
+
+void beyn_result_gsl_free(beyn_result_gsl_t *result);
+
+/// Beyn algorithm result structure (pure C array version).
+typedef struct beyn_result_t {
+	size_t neig; ///< Number of eigenvalues found.
+	size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec).
+	complex double *eigval;
+	complex double *eigval_err;
+	double *residuals;
+	complex double *eigvec;
+	double *ranktest_SV;
+
+	/// Private, we wrap it around beyn_result_gsl_t for now.
+	beyn_result_gsl_t *gsl;
+	
+} beyn_result_t;
+
+/// Converts beyn_result_gsl_t from beyn_result_t.
+/** After calling this, use beyn_result_free() on the returned pointer;
+ *  do NOT run beyn_result_gsl_free() anymore.
+ */
+beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g);
+
+void beyn_result_free(beyn_result_t *result);
+
+beyn_result_gsl_t *beyn_solve_gsl(
+		size_t m, ///< Dimension of the matrix \a M.
+		size_t l, ///< Number of columns of the random matrix  \f$ \hat V \f$ (larger than the expected number of solutions).
+		beyn_function_M_gsl_t M, ///< Function providing the matrix \f$ M(z) \f$.
+		beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, ///< Fuction providing the matrix \f$ M^{-1}(z) \hat V \f$ (optional).
+		void *params, ///< Parameter pointer passed to M() and M_inv_Vhat().
+		const beyn_contour_t *contour, ///< Integration contour.
+		double rank_tol, ///< (default: `1e-4`) TODO DOC.
+		double res_tol ///< (default: `0.0`) TODO DOC.
+	      );
+
+beyn_result_t *beyn_solve(
+		size_t m, ///< Dimension of the matrix \a M.
+		size_t l, ///< Number of columns of the random matrix  \f$ \hat V \f$ (larger than the expected number of solutions).
+		beyn_function_M_t M, ///< Function providing the matrix \f$ M(z) \f$.
+		beyn_function_M_inv_Vhat_t M_inv_Vhat, ///< Fuction providing the matrix \f$ M^{-1}(z) \hat V \f$ (optional).
+		void *params, ///< Parameter pointer passed to M() and M_inv_Vhat().
+		const beyn_contour_t *contour, ///< Integration contour.
+		double rank_tol, ///< (default: `1e-4`) TODO DOC.
+		double res_tol ///< (default: `0.0`) TODO DOC.
+	      );
+
+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/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd
index ba4df0c..0f808b1 100644
--- a/qpms/qpms_cdefs.pxd
+++ b/qpms/qpms_cdefs.pxd
@@ -589,3 +589,62 @@ cdef extern from "ewald.h":
             bint pgen_generates_shifted_points, cart3_t k, cart3_t particle_shift)
 
 
+cdef extern from "gsl/gsl_complex.h":
+    ctypedef struct gsl_complex:
+        double dat[2]
+
+cdef extern from "gsl/gsl_matrix.h":
+    ctypedef struct gsl_matrix_complex:
+        pass
+    ctypedef struct gsl_vector:
+        pass
+    ctypedef struct gsl_vector_complex:
+        pass
+
+cdef extern from "beyn.h":
+    ctypedef struct beyn_contour_t:
+        bint (*inside_test)(beyn_contour_t *, cdouble z)
+        pass
+    ctypedef struct beyn_result_gsl_t:
+        pass
+    ctypedef struct beyn_result_t:
+        size_t neig
+        size_t vlen
+        cdouble *eigval
+        cdouble *eigval_err
+        double *residuals
+        cdouble *eigvec
+        double *ranktest_SV
+        beyn_result_gsl_t *gsl
+    ctypedef enum beyn_contour_halfellipse_orientation:
+        BEYN_CONTOUR_HALFELLIPSE_RE_PLUS
+        BEYN_CONTOUR_HALFELLIPSE_IM_PLUS
+        BEYN_CONTOUR_HALFELLIPSE_RE_MINUS
+        BEYN_CONTOUR_HALFELLIPSE_IM_MINUS
+
+    ctypedef int (*beyn_function_M_gsl_t)(gsl_matrix_complex *target_M, cdouble z, void *params)
+    ctypedef int (*beyn_function_M_inv_Vhat_gsl_t)(gsl_matrix_complex *target, const gsl_matrix_complex *Vhat, cdouble z, void *params)
+    ctypedef int (*beyn_function_M_t)(cdouble *target_M, size_t m, cdouble z, void *params)
+    ctypedef int (*beyn_function_M_inv_Vhat_t)(cdouble *target, size_t m, size_t l, const cdouble *Vhat, cdouble z, void *params)
+
+    void beyn_result_gsl_free(beyn_result_gsl_t *result)
+    void beyn_result_free(beyn_result_t *result)
+
+    beyn_result_gsl_t *beyn_solve_gsl(size_t m, size_t l, beyn_function_M_gsl_t M,
+            beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, void *params, const beyn_contour_t *contour,
+            double rank_tol, double res_tol)
+
+    beyn_result_t *beyn_solve(size_t m, size_t l, beyn_function_M_t M,
+            beyn_function_M_inv_Vhat_t M_inv_Vhat, void *params, const beyn_contour_t *contour,
+            double rank_tol, double res_tol)
+
+    beyn_contour_t *beyn_contour_ellipse(cdouble centre, double halfax_re, double halfax_im, size_t npoints)
+    beyn_contour_t *beyn_contour_halfellipse(cdouble centre, double halfax_re, double halfax_im, size_t npoints,
+            beyn_contour_halfellipse_orientation ori)
+    beyn_contour_t *beyn_contour_kidney(cdouble centre, double halfax_re, double halfax_im, size_t npoints,
+            double rounding, beyn_contour_halfellipse_orientation ori)
+
+
+    cdouble gsl_comlpex_tostd(gsl_complex z)
+    gsl_complex gsl_complex_fromstd(cdouble z)
+
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 3d1ef0a..17ded6f 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -15,7 +15,68 @@ 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( kidneycontour kidneycontour.c )
+target_link_libraries( kidneycontour qpms gsl lapacke amos m )
+target_include_directories( kidneycontour PRIVATE .. )
+
+add_executable( tbeyn tbeyn.c )
+target_link_libraries( tbeyn qpms gsl lapacke amos m )
+target_include_directories( tbeyn PRIVATE .. )
+
+add_executable( tbeyn2 tbeyn2.c )
+target_link_libraries( tbeyn2 qpms gsl lapacke amos m )
+target_include_directories( tbeyn2 PRIVATE .. )
+
+add_executable( tbeyn3a tbeyn3.c )
+target_link_libraries( tbeyn3a qpms gsl lapacke amos m )
+target_include_directories( tbeyn3a PRIVATE .. )
+target_compile_definitions( tbeyn3a PRIVATE VARIANTA )
+
+add_executable( tbeyn3a_implus tbeyn3.c )
+target_link_libraries( tbeyn3a_implus qpms gsl lapacke amos m )
+target_include_directories( tbeyn3a_implus PRIVATE .. )
+target_compile_definitions( tbeyn3a_implus PRIVATE VARIANTA IMPLUS )
+
+add_executable( tbeyn3a_kidney tbeyn3.c )
+target_link_libraries( tbeyn3a_kidney qpms gsl lapacke amos m )
+target_include_directories( tbeyn3a_kidney PRIVATE .. )
+target_compile_definitions( tbeyn3a_kidney PRIVATE VARIANTA IMPLUS_KIDNEY )
+
+add_executable( tbeyn3b tbeyn3.c )
+target_link_libraries( tbeyn3b qpms gsl lapacke amos m )
+target_include_directories( tbeyn3b PRIVATE .. )
+target_compile_definitions( tbeyn3b PRIVATE VARIANTB RXSMALL )
+
+add_executable( tbeyn3bfail tbeyn3.c )
+target_link_libraries( tbeyn3bfail qpms gsl lapacke amos m )
+target_include_directories( tbeyn3bfail PRIVATE .. )
+target_compile_definitions( tbeyn3bfail PRIVATE VARIANTB )
+
+add_executable( tbeyn3c tbeyn3.c )
+target_link_libraries( tbeyn3c qpms gsl lapacke amos m )
+target_include_directories( tbeyn3c PRIVATE .. )
+target_compile_definitions( tbeyn3c PRIVATE VARIANTC RXSMALL )
+
+add_executable( tbeyn3d tbeyn3.c )
+target_link_libraries( tbeyn3d qpms gsl lapacke amos m )
+target_include_directories( tbeyn3d PRIVATE .. )
+target_compile_definitions( tbeyn3d PRIVATE VARIANTD RXSMALL )
+
+add_executable( tbeyn3e tbeyn3.c )
+target_link_libraries( tbeyn3e qpms gsl lapacke amos m )
+target_include_directories( tbeyn3e PRIVATE .. )
+target_compile_definitions( tbeyn3e PRIVATE VARIANTE RXSMALL )
+
+add_executable( tbeyn3f tbeyn3.c )
+target_link_libraries( tbeyn3f qpms gsl lapacke amos m )
+target_include_directories( tbeyn3f PRIVATE .. )
+target_compile_definitions( tbeyn3f PRIVATE VARIANTF )
+
+add_executable( tbeyn_gsl tbeyn_gsl.c )
+target_link_libraries( tbeyn_gsl qpms gsl lapacke amos m )
+target_include_directories( tbeyn_gsl 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/kidneycontour.c b/tests/kidneycontour.c
new file mode 100644
index 0000000..8f4e635
--- /dev/null
+++ b/tests/kidneycontour.c
@@ -0,0 +1,19 @@
+#include <qpms/beyn.h>
+#include <stdio.h>
+
+#define CPAIR(x) creal(x), cimag(x)
+
+int main(int argc, char **argv)
+{
+  double rRe = 2e3, rIm = 1.5e3, rounding = 0.2;
+  complex double centre = 1e3 * I;
+  size_t n = 100;
+
+  beyn_contour_t *c = beyn_contour_kidney(centre, rRe, rIm, rounding, n, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS);
+
+  for(size_t i = 0; i < n; ++i) 
+    printf("%g\t%g\t%g\t%g\n", CPAIR(c->z_dz[i][0]), CPAIR(c->z_dz[i][1]));
+
+  free(c);
+  return 0;
+}
diff --git a/tests/tbeyn.c b/tests/tbeyn.c
new file mode 100644
index 0000000..e2f52f4
--- /dev/null
+++ b/tests/tbeyn.c
@@ -0,0 +1,40 @@
+#include <qpms/beyn.h>
+#include <stdio.h>
+#include <string.h>
+
+// Matrix as in Beyn, section 4.11
+int M_function(complex double *target, const size_t m, const complex double z, void *no_params) {
+  complex double d =  2*m - 4*z / (6*m);
+  complex double od =  -((double)m) - z / (6*m);
+
+  memset(target, 0, m*m*sizeof(complex double));
+  for (int i = 0; i < m; ++i) {
+    target[i*m + i] = d;
+    if(i > 0) target[i*m + i-1] = od;
+    if(i < m - 1) target[i*m + i+1] = od;
+  }
+  target[m*(m-1) + m-1] = d/2 + z/(z-1);
+
+  return 0;
+}
+
+int main() {
+  complex double z0 = 150+2*I;
+  double Rx = 148, Ry = 148;
+  int L = 10, N = 50, dim = 400;
+  beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
+
+  beyn_result_t *result =
+    beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
+      contour, 1e-4, 1e-4);
+  printf("Found %zd eigenvalues:\n", result->neig);
+  for (size_t i = 0; i < result->neig; ++i) {
+    complex double eig = result->eigval[i];
+    printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig));
+  }
+  free(contour);
+  beyn_result_free(result);
+  return 0;
+}
+
+
diff --git a/tests/tbeyn2.c b/tests/tbeyn2.c
new file mode 100644
index 0000000..43b37d5
--- /dev/null
+++ b/tests/tbeyn2.c
@@ -0,0 +1,62 @@
+#include <qpms/beyn.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+
+static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); }
+// Normal distribution via Box-Muller transform
+static double randN(double sigma, double mu) {
+  double u1 = randU(0,1);
+  double u2 = randU(0,1);
+  return mu + sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2);
+}
+
+struct param49 {
+  double *T0;
+  double *T1;
+  double *T2;
+};
+
+// Matrix as in Beyn, example 4.9
+int M_function(complex double *target, const size_t m, const complex double z, void *params) {
+  struct param49 *p = params;
+
+  for(size_t i = 0; i < m*m; ++i)
+    target[i] = p->T0[i] + z*p->T1[i] + z*z*p->T2[i];
+  return 0;
+}
+
+int main(int argc, char **argv) {
+  complex double z0 = 0;
+  double Rx = .3, Ry = .3;
+  int L = 30, N = 150, dim = 60;
+  if (argc > 1) N = atoi(argv[1]);
+  beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
+  struct param49 p;
+  p.T0 = malloc(dim*dim*sizeof(double));
+  p.T1 = malloc(dim*dim*sizeof(double));
+  p.T2 = malloc(dim*dim*sizeof(double));
+  for(size_t i = 0; i < dim*dim; ++i) {
+    p.T0[i] = randN(1,0);
+    p.T1[i] = randN(1,0);
+    p.T2[i] = randN(1,0);
+  }
+
+  beyn_result_t *result =
+    beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, &p /*params*/,
+      contour, 1e-4, 1e-4);
+  printf("Found %zd eigenvalues:\n", result->neig);
+  for (size_t i = 0; i < result->neig; ++i) {
+    complex double eig = result->eigval[i];
+    printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig));
+  }
+  free(contour);
+  beyn_result_free(result);
+  free(p.T0);
+  free(p.T1);
+  free(p.T2);
+  return 0;
+}
+
+
diff --git a/tests/tbeyn3.c b/tests/tbeyn3.c
new file mode 100644
index 0000000..52ede9a
--- /dev/null
+++ b/tests/tbeyn3.c
@@ -0,0 +1,98 @@
+#define _GNU_SOURCE
+#include <qpms/beyn.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <fenv.h>
+
+static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); }
+// Normal distribution via Box-Muller transform
+static double randN(double sigma, double mu) {
+  double u1 = randU(0,1);
+  double u2 = randU(0,1);
+  return mu + sigma*sqrt(-2*log(u1))*cos(2.*M_PI*u2);
+}
+
+struct param {
+  double *T0;
+  double *T1;
+  double *T2;
+};
+
+int M_function(complex double *target, const size_t m, const complex double z, void *params) {
+  struct param *p = params;
+
+  for(size_t i = 0; i < m*m; ++i)
+    target[i] = p->T0[i] + z*p->T1[i] + cexp(
+#ifdef VARIANTB // Also note that this case requires pretty large contour point number (>~ 3000)
+        (1+3*I)
+#else // VARIANTA or VARIANTC
+        (1+1*I)
+#endif
+        *z*p->T2[i]) + 
+#ifdef VARIANTC // Essential singularity at zero
+      cexp(3/z);
+#elif defined VARIANTD // Essential singularity outside the contour
+      cexp(3/(z-1))
+#elif defined VARIANTE // High-order pole at zero
+      3/cpow(z,10);
+#elif defined VARIANTF // High-order pole at zero, higher order than dim
+      .0003/cpow(z,12);
+#else // double pole at zero
+      3/z/z
+#endif
+      ;
+  return 0;
+}
+
+int main(int argc, char **argv) {
+  feenableexcept(FE_INVALID | FE_OVERFLOW);
+  complex double z0 = 0+3e-1*I;
+#ifdef RXSMALL
+  double Rx = .1;
+#else 
+  double Rx = .3; // Variant B will fail in this case due to large number of eigenvalues (>30)
+#endif
+  double Ry = .25;
+#ifdef VARIANTF
+  int L = 10, N = 150, dim = 10;
+#else
+  int L = 30, N = 150, dim = 60;
+#endif
+  if (argc > 1) N = atoi(argv[1]);
+  if (argc > 2) L = atoi(argv[2]);
+#ifdef IMPLUS
+  beyn_contour_t *contour = beyn_contour_halfellipse(z0, Rx, Ry, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS);
+#elif defined IMPLUS_KIDNEY
+  beyn_contour_t *contour = beyn_contour_kidney(z0, Rx, Ry, 0.3, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS);
+#else
+  beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
+#endif
+  struct param p;
+  p.T0 = malloc(dim*dim*sizeof(double));
+  p.T1 = malloc(dim*dim*sizeof(double));
+  p.T2 = malloc(dim*dim*sizeof(double));
+  for(size_t i = 0; i < dim*dim; ++i) {
+    p.T0[i] = randN(1,0);
+    p.T1[i] = randN(1,0);
+    p.T2[i] = randN(1,0);
+  }
+
+  beyn_result_t *result =
+    beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, &p /*params*/,
+      contour, 1e-4, 1e-4);
+  printf("Found %zd eigenvalues:\n", result->neig);
+  for (size_t i = 0; i < result->neig; ++i) {
+    complex double eig = result->eigval[i];
+    printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig));
+  }
+  free(contour);
+  beyn_result_free(result);
+  free(p.T0);
+  free(p.T1);
+  free(p.T2);
+  return 0;
+}
+
+
diff --git a/tests/tbeyn_gsl.c b/tests/tbeyn_gsl.c
new file mode 100644
index 0000000..9f73b0c
--- /dev/null
+++ b/tests/tbeyn_gsl.c
@@ -0,0 +1,42 @@
+#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( -(double)m - z / (6*m) );
+
+  gsl_matrix_complex_set_zero(target);
+  for (int i = 0; i < m; ++i) {
+    gsl_matrix_complex_set(target, i, i, d);
+    if(i > 0) gsl_matrix_complex_set(target, i, i-1, od);
+    if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, od);
+  }
+  gsl_matrix_complex_set(target, m-1, m-1, gsl_complex_fromstd(gsl_complex_tostd(d)/2 + z/(z-1)));
+
+  return 0;
+}
+
+int main() {
+  complex double z0 = 150+2*I;
+  double Rx = 148, Ry = 148;
+  int L = 10, N = 50, dim = 400;
+  beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
+
+  beyn_result_gsl_t *result =
+    beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
+      contour, 1e-4, 1e-4);
+  printf("Found %zd eigenvalues:\n", result->neig);
+  for (size_t i = 0; i < result->neig; ++i) {
+    gsl_complex eig = gsl_vector_complex_get(result->eigval, i);
+    printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
+  }
+  free(contour);
+  beyn_result_gsl_free(result);
+  return 0;
+}
+
+