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};
clf
shg
set(gcf,'color','white')
for d = 1:2
    % Plot the region.
    vs = vertices{d};
    subplot(2,2,d)
    plot(vs(1,:),vs(2,:),'k.-');
    axis([-0.1 3.1 -0.1 3.1])
    axis square
    title(sprintf('drum%d',d));
end
figure(2)
% Generate a coarse finite difference grid.
ngrid = 5;
h = 1/ngrid;
[x,y] = meshgrid(0:h:3);
% Loop over the two regions.
for d = 1:2
    % Determine points inside and on the boundary.
    vs = vertices{d};
    [in,on] = inpolygon(x,y,vs(1,:),vs(2,:));
    in = xor(in,on);
    % Plot the region and the grid.
    subplot(2,2,d)
    plot(vs(1,:),vs(2,:),'k-',x(in),y(in),'b.',x(on),y(on),'k.')
    axis([-0.1 3.1 -0.1 3.1])
    axis square
    title(sprintf('drum%d',d));
    grid{d} = double(in);
end
figure(3)
for d = 1:2
    % Number the interior grid points.
    G = grid{d};
    p = find(G);
    G(p) = (1:length(p))';
    % Display the numbering.
    fprintf('grid%d =',d);
    spy(flipud(G))
    % Discrete Laplacian.
    A = delsq(G);
    % Spy plot
    subplot(2,2,d)
    markersize = 6;
    spy(A,markersize)
    title(sprintf('delsq(grid%d)',d));
end
%% Compute the eigenvalues

% How many eigenvalues?

   eignos = 20;

   % A finer grid.
   ngrid = 32;
   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);

   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

   format long
   lambda = flipud(lambda)