diff -r 46def32e6806 libinterp/corefcn/__eigs__.cc
--- a/libinterp/corefcn/__eigs__.cc	Mon Mar 29 07:54:26 2021 +0200
+++ b/libinterp/corefcn/__eigs__.cc	Tue Mar 30 21:12:30 2021 +0200
@@ -47,29 +47,31 @@
 
 #if defined (HAVE_ARPACK)
 
-// Global pointer for user defined function.
-static octave_value eigs_fcn;
-
-// Have we warned about imaginary values returned from user function?
-static bool warned_imaginary = false;
+struct EigsContext {
+  // Global pointer for user defined function.
+  octave_value eigs_fcn;
+
+  // Have we warned about imaginary values returned from user function?
+  bool warned_imaginary = false;
+};
 
 // Is this a recursive call?
 static int call_depth = 0;
 
 ColumnVector
-eigs_func (const ColumnVector& x, int& eigs_error)
+eigs_func (const ColumnVector& x, int& eigs_error, EigsContext& context)
 {
   ColumnVector retval;
   octave_value_list args;
   args(0) = x;
 
-  if (eigs_fcn.is_defined ())
+  if (context.eigs_fcn.is_defined ())
     {
       octave_value_list tmp;
 
       try
         {
-          tmp = octave::feval (eigs_fcn, args, 1);
+          tmp = octave::feval (context.eigs_fcn, args, 1);
         }
       catch (octave::execution_exception& ee)
         {
@@ -78,10 +80,10 @@
 
       if (tmp.length () && tmp(0).is_defined ())
         {
-          if (! warned_imaginary && tmp(0).iscomplex ())
+          if (! context.warned_imaginary && tmp(0).iscomplex ())
             {
               warning ("eigs: ignoring imaginary part returned from user-supplied function");
-              warned_imaginary = true;
+              context.warned_imaginary = true;
             }
 
           retval = tmp(0).xvector_value ("eigs: evaluation of user-supplied function failed");
@@ -97,19 +99,19 @@
 }
 
 ComplexColumnVector
-eigs_complex_func (const ComplexColumnVector& x, int& eigs_error)
+eigs_complex_func (const ComplexColumnVector& x, int& eigs_error, EigsContext& context)
 {
   ComplexColumnVector retval;
   octave_value_list args;
   args(0) = x;
 
-  if (eigs_fcn.is_defined ())
+  if (context.eigs_fcn.is_defined ())
     {
       octave_value_list tmp;
 
       try
         {
-          tmp = octave::feval (eigs_fcn, args, 1);
+          tmp = octave::feval (context.eigs_fcn, args, 1);
         }
       catch (octave::execution_exception& ee)
         {
@@ -196,9 +198,8 @@
   ColumnVector resid;
   ComplexColumnVector cresid;
   octave_idx_type info = 1;
-
-  warned_imaginary = false;
-
+  EigsContext context;
+
   octave::unwind_protect_var<int> restore_var (call_depth);
   call_depth++;
 
@@ -208,9 +209,9 @@
   if (args(0).is_function_handle () || args(0).is_inline_function ()
       || args(0).is_string ())
     {
-      eigs_fcn = octave::get_function_handle (interp, args(0), "x");
-
-      if (eigs_fcn.is_undefined ())
+      context.eigs_fcn = octave::get_function_handle (interp, args(0), "x");
+
+      if (context.eigs_fcn.is_undefined ())
         error ("eigs: unknown function");
 
       if (nargin < 2)
@@ -425,7 +426,11 @@
 
   octave_idx_type nconv;
   if (a_is_complex || b_is_complex)
-    {
+    {
+      EigsComplexFunc eigs_complex_func = [&context] (const ComplexColumnVector& x, int& eigs_error) {
+                                            return ::eigs_complex_func(x, eigs_error, context);
+                                          };
+
       ComplexMatrix eig_vec;
       ComplexColumnVector eig_val;
 
@@ -485,7 +490,11 @@
         }
     }
   else if (sigmai != 0.0)
-    {
+    {
+      EigsComplexFunc eigs_complex_func = [&context] (const ComplexColumnVector& x, int& eigs_error) {
+                                            return ::eigs_complex_func(x, eigs_error, context);
+                                          };
+
       // Promote real problem to a complex one.
       ComplexMatrix eig_vec;
       ComplexColumnVector eig_val;
@@ -534,6 +543,10 @@
     }
   else
     {
+      EigsFunc eigs_func = [&context] (const ColumnVector& x, int& eigs_error) {
+                             return ::eigs_func(x, eigs_error, context);
+                           };
+
       if (symmetric)
         {
           Matrix eig_vec;
diff -r 46def32e6806 liboctave/numeric/eigs-base.h
--- a/liboctave/numeric/eigs-base.h	Mon Mar 29 07:54:26 2021 +0200
+++ b/liboctave/numeric/eigs-base.h	Tue Mar 30 21:12:30 2021 +0200
@@ -28,6 +28,7 @@
 
 #include "octave-config.h"
 
+#include <functional>
 #include <iosfwd>
 #include <string>
 
@@ -38,10 +39,9 @@
 class Matrix;
 class ComplexMatrix;
 
-typedef ColumnVector (*EigsFunc) (const ColumnVector& x, int& eigs_error);
+typedef std::function<ColumnVector(const ColumnVector& x, int& eigs_error)> EigsFunc;
 
-typedef ComplexColumnVector (*EigsComplexFunc) (const ComplexColumnVector& x,
-                                                int& eigs_error);
+typedef std::function<ComplexColumnVector(const ComplexColumnVector& x, int& eigs_error)> EigsComplexFunc;
 
 template <typename M>
 OCTAVE_API octave_idx_type
