When we apply separation of variables to the Helmholtz equation of an elliptic membrane, we get Mathieu’s system. This two-parameter eigenvalue problem (2EP) is often used as a motivation for the introduction of multiparameter eigenvalue problems (MEPs), see, e.g., [41]. Yet, it was not until [16] that this approach was actually used to compute the eigenfrequencies of an elliptic membrane. In [16] it is shown that the two-parameter approach has certain advantages with respect to the accuracy as well as to the required computational time and can be used in practice to numerically evaluate a large number of eigenfrequencies.

Compared with [16], this paper presents advances in several directions. First, we consider several other very important problems. In addition to Mathieu’s system, discretization with spectral collocation in conjunction with numerical methods for the obtained algebraic MEP can be applied to Lamé’s system as well as other MEPs that arise by separation of variables. To the best of our knowledge, this technique

has never been applied to these problems until now. In particular, in Example 5 we apply this method to a three-parameter eigenvalue problem (3EP) and compute eigenmodes of a challenging ellipsoidal wave equation.

Second, while a Jacobi-Davidson style solver [19, 20] and standard implicitly restarted Arnoldi [37] were used in [16] to solve the 2EPs, we speed up the computations in this paper even more by exploiting a new fast Sylvester-Arnoldi algorithm for 2EPs from [29]; this algorithm is briefly described in Section 6. Third, we compare with, and improve on several available results in the literature in the numerical examples in Section 7. Finally, we present a MATLAB toolbox for solving MEPs coming from various equations.

The rest of this paper is organized as follows. Section 2 provides several examples of the technique of separation of variables leading to MEPs. A generic form of such MEP is treated in Section 3. In Section 4 we give an overview of the spectral collocation method, which is used to discretize a MEP into an algebraic MEP. This eigenvalue problem is presented in Section 5. In Section 6 we give an overview of available numerical methods for algebraic MEPs with an emphasis on the recent Sylvester-Arnoldi method from [29]. An important part of the paper are the numerical examples in Section 7. In several examples we demonstrate that spectral collocation combined with the Sylvester-Arnoldi method can compute several hundreds of the smallest eigenmodes very efficiently and accurately; we hereby improve various previous results. In the appendix we describe the main functions in a freely available MATLAB toolbox MultiParEig that contains the implementations of all algorithms and numerical examples from this paper.

Whenever the separation of variables is used to solve a boundary value problem related to a PDE, a system of ODEs is obtained in the first instance. Then, the boundary conditions of the problem at hand dictate boundary conditions for the unknowns of the systems of ODEs involved. Consequently, a MEP is now well defined.

In this section we give some examples of boundary value problems where separation of variables leads to MEPs. We do not attempt to describe all possible situations; for a good overview of all possible coordinate systems and related boundary value problems, see, e.g., [25, 30, 32, 46]. Additional examples together with numerical solutions can be found in Section 7.

An overview of the theory related to MEPs that are obtained by separation of variables is presented in [6]; see also [12, 36]. In the generic case separation of variables applied to a separable boundary value problem leads to a system of $k$ differential equations of the form

together with the appropriate boundary conditions, where $k=2$ or $k=3$. We are interested in a $k$-tuple $\left(\lambda_{1},\ldots,\lambda_{k}\right)$ and nontrivial functions $y_{1},\ldots,y_{k}$ such that equations (11) and the boundary conditions are satisfied. In such case we say that $\left(\lambda_{1},\ldots,\lambda_{k}\right)$ is an eigenvalue of (11). We remark that the numerical
methods from this paper are suitable for all problems of the form (11), regardless if they come from separation of variables or not.

If $x_{i}\in\left[a_{i},b_{i}\right]$, then generic boundary conditions are

where $\left(\alpha_{i},\beta_{i}\right)\neq(0,0)$ and $\left(\gamma_{i},\delta_{i}\right)\neq(0,0)$ for $i=1,\ldots,k$. In some cases boundary conditions can depend on eigenvalue parameters as well.

If the determinant $\delta\left(x_{1},\ldots,x_{k}\right):=\operatorname{det}\left(s_{ij}\left(x_{i}\right)\right)$ is definite for all $x_{i}\in\left[a_{i},b_{i}\right],i=1,\ldots,n$, then it is well-known that the problem (11) is solvable, see, e.g., [6, 41]. In addition, the eigenvalues can be ordered by the multi-indices $\left(n_{1},\ldots,n_{k}\right)$, which means that the corresponding eigenfunction $y_{i}\left(x_{i}\right)$ has exactly $n_{i}$ zeros on $\left(a_{i},b_{i}\right)$ for $i=1,\ldots,k$.

There exist some numerical methods for MEPs of the form (11). For problems that have additional definiteness properties, a continuation method is proposed in [2]; see also [23]. Once the problem (11) is discretized into an algebraic MEP, all methods from Section 6 can be applied.

Most of the numerical methods are limited to 2EPs, i.e., $k=2$. A generalization of the shooting method is presented in [7] and a two-dimensional bisection is presented in [22]. A method using the Prüfer transformation that can compute an eigenvalue with a prescribed multi-index is presented in [1] for 2EPs and generalized to 3EPs, i.e., $k=3$, in [26]. For some numerical methods for 3EPs see, e.g., [27] and the references therein.

In this section we briefly introduce the spectral collocation method, which requires smaller matrices and returns more accurate results when compared with finite differences and finite elements methods. For more details, see, e.g., [10, 13, 15]. The spectral collocation method is based on approximating the solution $y(x)$ with a finite linear combination of a chosen set of orthogonal function as

The coefficients $\alpha_{1},\ldots,\alpha_{N}$ are chosen so that $y\left(x_{i}\right)=y_{N-1}\left(x_{i}\right)$ for $i=1,\ldots,N$, where $x_{1},\ldots,x_{N}$ are distinct collocation nodes. In this sense, $y_{N-1}(x)$ is an interpolating function for $y(x)$ on the set of collocation nodes. The large positive integer $N$ is the order of approximation and is also called the cutoff parameter. To accommodate with physics, $N$ can be viewed as the number of degrees of freedom (of modes) of a system.

Depending on the interval and other properties of the problem, various combinations of basis functions and collocation nodes exist. If the interval is finite, then a common choice is the set of Chebyshev polynomials as basis functions and collocation points are the Chebyshev nodes of the second kind. The method is called Chebyshev collocation ( ChC ) and the nodes on the interval [ $-1,1$ ] are defined by

for $j=1,\ldots,N$ (see, e.g., [42, Eq. (13)]).
As we have to discretize each of the equations of (11) independently, it is sufficient to consider one equation. Suppose that we have a differential equation

on the interval $[a,b]$. The ChC discretization of the above equation is

where $D_{N}$ and $D_{N}^{2}$ are first and second order differentiation matrices in the Chebyshev nodes of the second kind and $f=\left[y\left(x_{1}\right),\ldots,y\left(x_{N}\right)\right]^{T}$ is the vector of values in the collocation nodes

for $j=1,\ldots,N$, which are translated from $[-1,1]$ to $[a,b]$. The matrices $P,Q,R,S$, and $T$ are diagonal with diagonal elements $p\left(x_{j}\right),q\left(x_{j}\right),r\left(x_{j}\right),s\left(x_{j}\right)$, and $t\left(x_{j}\right)$, respectively

We use the seminal paper [43] to obtain the entries of the differentiation matrices $D_{N}$ and $D_{N}^{2}$; see also [11] for Maple code and [40] for MATLAB code. The matrices $D_{N}$ and $D_{N}^{2}$ are dense, non-symmetric, and have much higher condition numbers as the corresponding symmetric tridiagonal finite difference matrices of the same size; see, e.g., [15, Fig. 3.13]. However, since discretization with the ChC requires much smaller matrices than with finite differences (see, e.g., [16, Ex. 5]), $N$ is usually small enough so that the ill-conditioning is not an issue.

To impose the boundary conditions, we apply the simple and general strategy from [21]. Another option is used if boundary conditions depend on eigenvalue parameters. In such case we replace the first and the last equation of (14) with the equations from the boundary conditions. For instance, suppose that the boundary condition for (13) in point $a$ is

We write the above as

where $D_{N}(1,:)$ is the first row of $D_{N}$, and replace the first row of (14) with the above equation. We proceed similarly in point $b$.

In many cases the boundary conditions are behavioral, which means that they are satisfied implicitly by choosing basis functions that have the required property or behavior. For instance, suppose that the differential equation (13) has singularity at the endpoint $a$ such that $p(a)=0$ and the boundary condition is that the solution should be bounded at $a$. As the Chebyshev polynomials (translated from $[-1,1]$ to $[a,b]$ ) are analytic at the endpoints, their sum satisfies the boundedness condition at $x=a$ and no additional explicit conditions are needed. Since collocation points in the ChC include the endpoints, the first equation of (14), which corresponds to the collocation point $x=a$, imposes the boundary condition at $x=a$ obtained from (13) by taking the limit $x\rightarrow a$. For more details, see [10, Section 6.3].

In case of the half line [ $0,\infty$ ) we use the Laguerre collocation (LC) (see, e.g., [15]). The basis functions are products of Laguerre polynomials and the weight function $e^{-\frac{1}{2}bx}$, where $b>0$ is the scaling parameter. The collocation nodes are $x_{1}=0<x_{2}/b<\cdots<x_{N}/b$, where $x_{2},\ldots,x_{N}$ are the roots of the Laguerre polynomial $L_{N-1}$ of degree $N-1$. When we use the LC, the boundary condition at infinity, which has to be $\lim_{x\rightarrow\infty}y(x)=0$, is automatically satisfied. The LC differentiation is exact for functions of the form $y(x)=e^{-\frac{1}{2}bx}p(x)$, where $p(x)$ is a polynomial of degree $\leq N-1$. The scaling parameter $b$ has to be carefully chosen according to how fast the solution decays to zero as $x\rightarrow\infty$.

When we discretize the MEP (11) we obtain an algebraic MEP, which has the form

where $A_{ij}$ is an $n_{i}\times n_{i}$ complex matrix for $i=1,\ldots,k$ and $j=0,\ldots,k$. A $k$-tuple ( $\lambda_{1},\ldots,\lambda_{k}$ ) is an eigenvalue if it satisfies (15) for nonzero vectors $x_{1}\in\mathbb{C}^{n_{1}},\ldots,x_{k}\in\mathbb{C}^{n_{k}}$. The corresponding eigenvector is the tensor product $x_{1}\otimes\cdots\otimes x_{k}$. Introducing the so-called $k\times k$ operator determinants

for $i=1,\ldots,k$, where the Kronecker product $\otimes$ is used instead of multiplication, we obtain matrices $\Delta_{0},\ldots,\Delta_{k}$ of size $n_{1}\cdots n_{k}\times n_{1}\cdots n_{k}$. If $\Delta_{0}$ is nonsingular, then the matrices $\Delta_{0}^{-1}\Delta_{1},\ldots,\Delta_{0}^{-1}\Delta_{k}$ commute and (15) is equivalent (for details, see, e.g., [5]) to a system of generalized eigenvalue problems

for decomposable tensors $z=x_{1}\otimes\cdots\otimes x_{k}$. This relation enables one to use standard numerical methods for the computation of eigenvalues if the $\Delta$-matrices are not too large. However, when we use spectral methods to discretize (11), then usually even for $k=2$ the $\Delta$-matrices are so large that it is not efficient or even not feasible to compute all the eigenvalues. The available numerical methods are presented in the next section.

Let us mention that if we discretize a self-adjoint MEP (11) obtained via the separation of variables using the spectral collocation method from Section 4, then the algebraic MEP is solvable, although it is not self-adjoint (since matrices $D_{N}$ and $D_{N}^{2}$ are not symmetric). Namely, the matrices $A_{ij}$ for $i,j=1,\ldots,k$ are all diagonal and if $\delta\left(x_{1},\ldots,x_{k}\right):=\operatorname{det}\left(s_{ij}\left(x_{i}\right)\right)$ is definite for all $x_{i}\in\left[a_{i},b_{i}\right],i=1,\ldots,n$, then the corresponding operator determinant $\Delta_{0}$ is nonsingular.

Example 1. For future reference, we give (15) and the corresponding operator determinants for the case $k=2$. An algebraic 2EP has the form

The eigenvalue is a pair ( $\lambda,\mu$ ) which satisfies (17) for nonzero vectors $x_{1}$ and $x_{2}$. The corresponding $2\times 2$ operator determinants of size $n_{1}n_{2}\times n_{1}n_{2}$ are

If $\Delta_{0}$ is nonsingular, then $\Delta_{0}^{-1}\Delta_{1}$ and $\Delta_{0}^{-1}\Delta_{2}$ commute and (17) is equivalent to a coupled pair of generalized eigenvalue problems

for decomposable tensors $z=x_{1}\otimes x_{2}$. This will be used in all our numerical examples.

If the size $n_{1}\cdots n_{k}$ of the $\Delta$-matrices in (16) is small enough, say $n_{1}\cdots n_{k}\leq 1000$ on a typical laptop, we can efficiently compute all eigenvalues of a MEP (15) from the related system of generalized eigenvalue problems (16). This numerical algorithm, which is based on the QZ algorithm, is given in Algorithm 1. It was presented in [19] for $k=2$, and can be generalized in a straightforward way to three or more parameters. The MATLAB toolbox MultiParEig that we describe in the appendix contains implementations of Algorithm 1 and its generalization for $k=3$.

Algorithm 1 An algorithm for the nonsingular two-parameter eigenvalue problem (17).

1. 1.

   Compute a generalized Schur decomposition $Q^{H}\Delta_{0}Z=R$ and $Q^{H}\Delta_{1}Z=S$, such that $Q$ and $Z$ are unitary, $R$ and $S$ are upper triangular, and multiple values of $\lambda_{i}:=s_{ii}/r_{ii}$ are clustered along the diagonal of $R^{-1}S$. As a result, $R$ and $S$ are partitioned as

   Report issue for preceding element

where the size of $R_{ii}$ and $S_{ii}$ is $m_{i}$ and $m_{1}+\cdots+m_{p}=n_{1}n_{2}$.
2. Compute diagonal blocks $T_{11},\ldots,T_{pp}$ of $T=Q^{H}\Delta_{2}Z$.
3. Compute the eigenvalues $\mu_{i1},\ldots,\mu_{im_{i}}$ of $T_{ii}w=\mu R_{ii}w$ for $i=1,\ldots,p$.
4. The eigenvalues of (17) are $\left(\lambda_{1},\mu_{11}\right),\ldots,\left(\lambda_{1},\mu_{1m_{1}}\right);\ldots;\left(\lambda_{p},\mu_{p1}\right),\ldots,\left(\lambda_{p},\mu_{pm_{p}}\right)$.
5. For each eigenvalue $(\lambda,\mu)$ compute the eigenvector $x_{1}\otimes x_{2}$ by inverse iteration on
$A_{i}-\lambda B_{i}-\mu C_{i}$ for $i=1,2$.

If $n_{1}\cdots n_{k}$ is too large then it is not efficient or even not feasible to compute all the eigenvalues. Instead, we can compute a subset of eigenvalues by an iterative method. For $k=2$, if we are looking for eigenvalues $(\lambda,\mu)$ close to a given target $(\sigma,\tau)$, then a method of choice is the harmonic version [20] of the Jacobi-Davidson method [19], which can also be easily generalized to three or more parameters. For an overview of other numerical methods, see, e.g., [19] and references therein.

In the remainder we will focus on numerical methods for 2EPs. Recently, a set of algorithms was presented in [29] that are very suitable for 2EPs related to separable boundary value problems. In many applications, for instance when we consider the Helmholtz equation $\nabla^{2}u+\omega^{2}u=0$, only one of the eigenvalue parameters $\lambda$ and $\mu$ is related to $\omega$ and thus relevant. Without loss of generality we can assume that this parameter is $\mu$. In a typical situation we are interested in a number $m$ of the lowest eigenfrequencies $\omega$ and the corresponding eigenmodes. Usually, the problem can be reformulated in such way that we have to find the $m$ eigenvalues $(\lambda,\mu)$ with the smallest $|\mu|$.

While a problem (11) usually has infinitely many eigenvalues, each discretization (17) has only finitely many eigenvalues. We may expect that some of them are good approximations to the eigenvalues of (11) while the remaining ones contain no useful information. For one-parameter problems, a heuristic suggested in [10, Rule-of-Thumb 8, p. 132] states:

In solving a linear eigenvalue problem by a spectral method using $N+1$ terms in the truncated spectral series, the lowest $\frac{1}{2}N$ eigenvalues are usually accurate to within a few percent while the larger $\frac{1}{2}N$ numerical eigenvalues differ from those of differential equation by such large amounts as to be useless.

Much time and energy has been spent along the last decades in the invention of slight modifications to standard spectral methods in order to solve the "spurious eigenvalue" difficulty (see, e.q., [18] and [10, Sections 7.5 and 7.6]). On one hand there are physically spurious eigenvalues that are numericallycomputed ones which are in error because of misapplication of boundary conditions or some other misrepresentation of the physics. The original continuous MEPs that we consider are self-adjoint, physically well established, and consequently the appearance of such spurious eigenvalues is highly unlikely. On the other hand we can always expect some numerically spurious eigenvalues that are poor approximation to exact eigenvalues because the mode is oscillating too rapidly to be resolved by $N$ degrees of freedom. Such numerically spurious eigenvalue can be computed accurately by using sufficiently large $N$. This explains the fairly high accuracy in computing the "smallest" eigenvalue for a self-adjoint problem with target oriented algorithms.

As a reliable test it is suggested in [10] to repeat the calculations with different $N$ and compare the results. It is also noted that the number of good eigenvalues may be smaller if the modes have boundary layers, critical levels, or other areas of very rapid change, or when the interval is unbounded. We will not discuss the convergence here, but, as we increase the sizes of the matrices in (17), more and more eigenvalues of (17) should converge to exact eigenvalues of (11) and the eigenvalues with the smallest value of $|\mu|$ typically converge among first.

The sizes $n_{1}$ and $n_{2}$ depend on the number $m$ of the wanted lowest eigenfrequencies and on the required accuracy. Usually $m\ll n_{1}n_{2}$ and if $n_{1}n_{2}$ is small, we can compute all values $\mu$ from the generalized eigenvalue problem $\Delta_{2}z=\mu\Delta_{0}z$. If $n_{1}n_{2}$ is too large for Algorithm 1, then we can apply the implicitly restarted Arnoldi method [37], Krylov-Schur method [38], or any other iterative method based on Krylov subspaces, to the matrix $\Gamma_{2}^{-1}:=\Delta_{2}^{-1}\Delta_{0}$. To do this, we must be able to efficiently multiply by $\Gamma_{2}^{-1}$, i.e., to solve a system

for a given vector $z$. In principle, the complexity of the above step is $\mathscr{O}\left(n_{1}^{3}n_{2}^{3}\right)$ as the matrices $\Delta_{0},\Delta_{1}$, and $\Delta_{2}$ are of size $n_{1}n_{2}\times n_{1}n_{2}$. But, by formulating (20) as a Sylvester equation, we can solve it in $\mathscr{O}\left(n_{1}^{3}+n_{2}^{3}\right)$. The key to the big reduction is the well-known equality

which enables us to write (20), i.e., $\left(B_{1}\otimes A_{2}-A_{1}\otimes B_{2}\right)w=\left(B_{1}\otimes C_{2}-C_{1}\otimes B_{2}\right)z$, as

where $z=\operatorname{vec}(Z)$ and $w=\operatorname{vec}(W)$. We get the Sylvester equation

which can be solved in $\mathscr{O}\left(n_{1}^{3}+n_{2}^{3}\right)$ by the Bartels-Stewart algorithm [8]. The overall algorithm from [29], which is included in MATLAB toolbox MultiParEig, is given in Algorithm 2.

Here we give some additional explanation about Algorithm 2, for more details see [29]. The algorithm requires that $A_{1}$ and $A_{2}$ are both nonsingular. If this is not the case, then it follows from [29, Lemma 3.1] that there exists $\theta$ such that $A_{1}-\theta B_{1}$ and $A_{2}-\theta B_{2}$ are both nonsingular. When we find such $\theta$, which is not difficult as there are only finitely many inappropriate values, we apply the shift $\lambda=\widetilde{\lambda}+\theta$ and use Algorithm 2 on the shifted problem

1. 1.

   Compute Schur decompositions $A_{2}^{-1}B_{2}=U_{1}R_{1}U_{1}^{H}$ and $B_{1}^{T}A_{1}^{-T}=U_{2}R_{2}U_{2}^{H}$.

   Report issue for preceding element
2. 2.

   Apply implicitly restarted Arnoldi or Krylov-Schur method on $\Gamma_{2}^{-1}:=\Delta_{2}^{-1}\Delta_{0}$, where in each step the product $w=\Gamma_{2}^{-1}z$ is computed as:
   (a) set matrix $Z$ such that $z=\operatorname{vec}(Z)$
   (b) solve Sylvester equation $R_{1}W-WR_{2}=-U_{1}^{H}A_{2}^{-1}ZA_{1}^{-T}U_{2}$
   (c) $w=\operatorname{vec}(W)$

   Report issue for preceding element

Once we have an eigenpair ( $\mu,z$ ) of $\Delta_{2}z=\mu\Delta_{0}z$, it is easy to compute the eigenvalue ( $\lambda,\mu$ ) and the vectors $x_{1}$ and $x_{2}$ of (17). The eigenvector $z=x_{1}\otimes x_{2}$ is decomposable and if $z=\operatorname{vec}(Z)$, then $Z$ has rank one and $Z=x_{2}x_{1}^{T}$. Thus, we can extract $x_{1}$ and $x_{2}$ from the singular value decomposition of $Z$. From $x_{1}$ and $x_{2}$ the eigenvalue parameter $\lambda$ can be computed from the tensor Rayleigh quotient [34] efficiently as

In the case of three or more parameters, we are not aware of a method that reduces the complexity of solving the linear system with the $\Delta$-matrices in a similar way as in Algorithm 2; this is left for future research. We can apply any iterative method based on Krylov subspaces on the $\Delta$-matrices, but since we use full vectors of size $n_{1}\cdots n_{k}$, we cannot use as many collocation points as for 2EPs.

In this section we give some numerical examples that show the strengths of the proposed approach and explain some further details. From the literature we take several boundary value problems that lead to some 2EPs and one 3EP, and solve them numerically with the methods proposed in the paper using MATLAB toolbox MultiParEig [35]. In many cases our numerical results are more accurate than the previously reported in the literature and we also give more eigenmodes. Fast computational times show that it is possible to numerically solve many separable boundary value problems efficiently as MEPs.

The following numerical examples were obtained on 64-bit Windows version of Matlab R2012b running on $\operatorname{Intel}(\mathrm{R})\operatorname{Core}(\mathrm{TM})$ i5-4670 3.40 Ghz processor and 16 GB of RAM.

Example 1. In the first example, which is based on [31], we consider Lamé’s system (see Subsection 2.2), related to the problem of computing the strength of a charge singularity of a flat plate for corner angle $0<\chi<2\pi$ :

where $r\geq 0,0\leq\theta\leq 2\pi,0\leq\phi\leq\pi,0\leq k,k^{\prime}\leq 1$, and $k^{2}+k^{\prime 2}=1$. Since the problem is located in the half-space $z\geq 0$, the range of $\theta$ can be reduced to $[0,\pi]$. For the solution $N(\theta)$ it either holds $N(\boldsymbol{\pi}-\boldsymbol{\theta})=N(\boldsymbol{\theta})$ so that $N^{\prime}(\boldsymbol{\pi}/2)=0$, or $N(\boldsymbol{\pi}-\boldsymbol{\theta})=-N(\boldsymbol{\theta})$ and hence $N(\boldsymbol{\pi}/2)=0$. Therefore, it suffices to consider $\theta\in\left[0,\frac{\pi}{2}\right]$. The relevant boundary conditions for our case are $L(0)=L^{\prime}(\pi)=0$ and $N^{\prime}(0)=N^{\prime}(\pi/2)=0$ if $0<\chi<\pi$ or $N(0)=N^{\prime}(\pi/2)=0$ if $\pi<\chi<2\pi$. The goal is to compute the smallest $\rho>0$ for $k=\sin(|\pi-\chi|/2)$. This problem is solved numerically as a 2EP in [29] using finite
differences, where quite large matrices (of size $40000\times 40000$ ) are needed for accurate results. Using the ChC discretization more accurate results can be obtained with much smaller matrices.

We set $\lambda=\delta,\mu=\rho(\rho+1)$ and discretize (22) and (23) using the ChC. In the discretized problem (17), $A_{1}$ and $A_{2}$ are full, $C_{1}$ and $C_{2}$ are diagonal, and $-B_{1}$ and $B_{2}$ are identity matrices. Due to the boundary conditions, matrices $A_{2},C_{1}$, and $C_{2}$ are singular. We fix this by shifting $\lambda$ so that we can apply the implicitly restarted Arnoldi version of Algorithm 2. We take the rather arbitrary shift $\lambda=\widetilde{\lambda}-10$ and then numerically solve the shifted problem (21) for $\theta=-10$. Numerical experiments show that the use of 60 collocation points for each equation is sufficient to compute the ten smallest values $\rho$ to at least 8 accurate digits. All values, when rounded, agree to the values in [31, Tabs. 1 and 2], which are given with 4 or 5 digits. The same problem is solved in [7] with the shooting method and in [22] with the two-dimensional bisection method, they both report the corresponding values $\mu$. In Table 1 we give values for some corner angles $\chi$ for comparison. The remaining values can be computed using the function demo_lame in toolbox MultiParEig. Let us note that the computation of all values $\rho$ and $\mu$ (for 31 different values of $\chi$ ) takes just 0.6 seconds.