function [Ex, Ey, K2] = ComputeTMModes(NumModes, Nx, Ny, Dx, Dy)
% [Ex, Ey, K] = ComputeTMModes(NumModes, Nx, Ny, Dx, Dy)
%
% Solves the 2D eigenvalueproblem with Ez computed in the corners of
% each cell. Homogeneous dirichlet BC.
% size(Ex) = [Nx, Ny-1, NumModes] 
% size(Ey) = [Nx-1, Ny, NumModes] 
% size(K)  = [NumModes] 

  % Number of unknowns
  nx = Nx - 1;
  ny = Ny - 1;
  N = nx * ny;
  A = repmat(-2*(Dx^-2 + Dy^-2), N, 1);
  % X neighbours
  Bu = repmat(Dx^-2, N, 1);
  Bl = Bu;
  Bu( 1:nx:N) = 0;
  Bl(nx:nx:N) = 0;
  % Y neighbours
  C = repmat(Dy^-2, N, 1);
  % The total matrix
  M = spdiags([C Bl A Bu C], [-nx -1 0 1 nx], N, N);
  M = full(M);
  [V, K] = eig(-M);
  [K,I] = sort(diag(K));

  % Init outdata vectors
  Ex = zeros(Nx * (Ny - 1), NumModes);
  Ey = zeros((Nx - 1) * Ny, NumModes);
  K2 = zeros(1, NumModes);
  
  for k = 1:NumModes
    Ez      = reshape(V(:,I(k)), nx, ny);
    tmp     = diff([zeros(1,ny); Ez; zeros(1,ny)], 1, 1) / Dx;
    Ex(:,k) = tmp(:);
    tmp     = diff([zeros(nx,1), Ez, zeros(nx,1)], 1, 2) / Dy;
    Ey(:,k) = tmp(:);
    K2(k)     = K(I(k));
  end
