How to compute the first eigenvalue of Laplace operator in an ellipse?

I like this question! I didn't have the stamina to look into the issue too precisely, but here's what I got so far from Abramowitz and Stegun. As you said, the eigenfunction will be described by Mathieu functions. More precisely, on page 722 they study the "wave equation for the elliptic cylinder" (I would call that the Helmholtz equation...) using separation of variables. For your case you want to take $\varphi = $ constant because you have only two variables. This means that your eigenfunction is given by $f(u)g(v)$ where $f$ and $g$ solve the modified Mathieu and Mathieu equation respectively as in 20.1.6. In the equations you have a parameter $q$ which depends on the ellipse (I think for you $q = \lambda_1(c) \cdot c^2/4$) and you have a parameter $a$ which comes from separation of variables and which we must determine along with $\lambda_1(c)$. You can also find these equations on the wikipedia page for Mathieu functions.

Here $u$ is the "radial" variable and $v$ is the "angular" variable (the picture in the wikipedia page for elliptical coordinates which you linked to above depicts these nicely). This means that $g$ must be $2\pi$ periodic to define a function on the ellipse, and $f$ must vanish at $1$ to satisfy the Dirichlet boundary condition (to simplify this part I've renormalized your area so the ellipse passes through $(x,y) = (c\cosh 1,0)$). Also, both must be positive inside the ellipse because the ground state is always positive. Over the next few pages you have a study telling you which $a$ lead to periodic solutions. The only one which is never vanishing is $a_0(q)$ and this gives a solution $g(v) = \textrm{ce}_0(v)$ which is graphed in figure 20.2.

Now we've identified the angular part of the eigenfunction, namely $g(v) = \textrm{ce}_0(v)$, and the parameter $a = a_0(q) \approx - q^2/2$ (for $q$ small -- this is equivalent to $c$ being small) which has an expansion given in 20.2.25. To fix the eigenvalue $\lambda_1(c)$ we plug this into the equation for $f(u)$, and we see what is the smallest positive value which gives us a solution which has $f'(0) = 0$ (so the solution is smooth) and has $f(1) = 0$.

If my calculation is right, $f(\pm iv)$ solves the same equation as $g(v)$, so we may take (this part is not precise! -- you will need to make some kind of approximation argument here) $f(u) = \textrm{ce}_0(iu)$. Then the expansion 20.2.27 says that if $q$ is small then $f(u) = \textrm{ce}_0(iu) \approx 1 - \frac q 2 \cosh(2u)$ which has $f'(0) = 0$. If $f(1) = 0$ then $q = 2/\cosh2$, and we get $\lambda_1(c) = \frac 8 {c^2 \cosh(2)}$. To compare with your answer for a circle, we must renormalize the area. This ellipse has area $A =\pi c^2 \sinh 1 \cosh 1$, giving $$ \lambda_1(c) \approx \frac {8 \pi \sinh 1 \cosh 1}{A \cosh 2} \approx 3.9 \frac \pi A$$

This is not the right answer (should have gotten 5.8 rather than 3.9), and I think the issue here comes from the fact that the approximation $\textrm{ce}_0(iu) \approx 1 - \frac q 2 \cosh(2u)$ is not good near $u=1$, even if $q$ is very small. What you really need is info on the zeros of Mathieu functions, which I'm too tired to look up... do post again if you find anything!