 % Vertices
   drum1 = [0 0 2 2 3 2 1 1 0
            0 1 3 2 2 1 1 0 0];
   drum2 = [1 0 0 2 2 3 2 1 1
            0 1 2 2 3 2 1 1 0];
   vertices = {drum1,drum2};

   % Number of eigenvalues

   eignos = 20;

   % Grid size
   ngrid = 32;

   % Compute the eigenvalues and eigenfunctions

   h = 1/ngrid;
   [x,y] = meshgrid(0:h:3);

   inpoints = (7*ngrid-2)*(ngrid-1)/2;
   lambda = zeros(eignos,2);
   V = zeros(inpoints,eignos,2);

   % Loop over the two drums

   for d = 1:2
      vs = vertices{d};
      [in,on] = inpolygon(x,y,vs(1,:),vs(2,:));
      in = xor(in,on);

      % Number the interior grid points.

      G = double(in);
      p = find(G);
      G(p) = (1:length(p))';
      grid{d} = G;

      % The discrete Laplacian

      A = delsq(G)/h^2;

      % Sparse matrix eigenvalues and vectors.

      [V(:,:,d),E] = eigs(A,eignos,0);
      lambda(:,d) = diag(E);
   end

   % Plot the eigenfunctions.

   for d = 1:2
      for k = 1:eignos
         figure(ceil(k/2))
         set(gcf,'color','white')
         subplot(2,2,2*mod(k-1,2)+d)

         % Insert the k-th eigenvector in the grid interior.

         G = grid{d};
         p = find(G);
         u = zeros(size(G));
         u(p) = V(:,eignos+1-k,d);

         % Make first eigenvector positive so its color matches the others.

         if k == 1
            u = -u;
         end

         % Insert NaN's to make the exterior disappear.

         vs = vertices{d};
         [in,on] = inpolygon(x,y,vs(1,:),vs(2,:));
         u(~in) = NaN;

         % A filled contour plot with a line on the boundary.

         s = max(abs(u(:)));
         contourf(x,y,u,s*(-1:1/4:1))
         line(vs(1,:),vs(2,:),'color','black','linewidth',2)
         title(num2str(k))
         axis([-0.1 3.1 -0.1 3.1])
         axis square off
      end
   end

   
   