function [Ex, Ey, K2] = ComputeTEModes(NumModes, Nx, Ny, Dx, Dy)
% [Ex, Ey, K] = ComputeTEModes(NumModes, Nx, Ny, Dx, Dy)
%
% Solves the 2D eigenvalueproblem with Hz computed in the mid-point of
% each cell. Homogeneous Neumann BC.
% size(Ex) = [Nx, Ny-1, NumModes] 
% size(Ey) = [Nx-1, Ny, NumModes] 
% size(K)  = [NumModes] 

  % Number of unknowns
  N = Nx * Ny;
  % Main diagonal
  A = repmat(-2*(Dx^-2 + Dy^-2), N, 1);
  I = [1:Nx N-Nx+1:N];  % Compensate for BC
  A(I) = A(I) + Dy^-2; 
  I = [1:Nx:N Nx:Nx:N]; 
  A(I) = A(I) + Dx^-2; 
  % 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
    Hz      =  reshape(V(:,I(k+1)), Nx, Ny);
    tmp     =  diff(Hz, 1, 2) / Dy;
    Ex(:,k) =  tmp(:);
    tmp     = -diff(Hz, 1, 1) / Dx;
    Ey(:,k) =  tmp(:); 
    K2(k) = K(I(k+1));
  end
  


