% Eigenmodes of Laplacian on the unit square
clear all; clc;

% define the spatial domain and temporal domain
a = 0;
b = 1;

% spatial and temporal grid parameters
n = 101;
h = (b-a)/(n-1);

% generate the centered difference differentiation matrix for the Laplacian
e = 1/h/h*ones(n^2,1);
D2 = spdiags([e,e,-4*e,e,e],[-n,-1:1,n],n^2,n^2);

% define an index function for this grid
index_function = @(i,j)((j-1)*n+i);

% get the list of Dirichlet indices
non_dirichlet_indices = [];
dirichlet_indices = [];
for i=1:n
    for j=1:n
        if (i~=1) && (i~=n) && (j~=1) && (j~=n)
            non_dirichlet_indices = ...
                [non_dirichlet_indices,index_function(i,j)];
        else
            dirichlet_indices = ...
                [dirichlet_indices,index_function(i,j)];
        end
    end
end

% strip out the Dirichlet rows and columns
A = -D2(non_dirichlet_indices,non_dirichlet_indices);

n_eigs = 6;
[V,D] = eigs(A,n_eigs,'SM');

% get the meshgrid, for plotting
x = linspace(a,b,n);
y = linspace(a,b,n);

% initialize VV
xx = zeros(n,n);
yy = zeros(n,n);
VV = 0*xx;

eigen_index = n_eigs;

% construct VV so that we can plot the eigenmode
% note that we have to pad V with zeros to show
% the Dirichlet boundary conditions
count = 1;
for i=1:n
    for j=1:n
        xx(i,j) = x(i);
        yy(i,j) = y(j);
        
        index = index_function(i,j);
        if isempty(find(dirichlet_indices == index))
            VV(i,j) = V(count,eigen_index);
            count = count + 1;
        else
            VV(i,j) = 0;
        end
    end
end

% plot the solution
contour(xx,yy,VV,30,'linewidth',2);
title(['$\lambda = $',num2str(D(eigen_index,eigen_index))], ...
    'interpreter','latex','fontsize',18)
xlabel('$x$','interpreter','latex','fontsize',18)
ylabel('$y$','interpreter','latex','fontsize',18)