diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ac7a18e..8185d79 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -23,6 +23,41 @@ 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( 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 .. )
diff --git a/tests/tbeyn3.c b/tests/tbeyn3.c
new file mode 100644
index 0000000..0c983b6
--- /dev/null
+++ b/tests/tbeyn3.c
@@ -0,0 +1,89 @@
+#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 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) {
+  complex double z0 = 0;
+#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 = .3;
+#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]);
+  beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
+  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;
+}
+
+