% tested with GNU Octave, version 3.6.2-5+deb7u1 x86_64-pc-linux-gnu

% number of modes to compute
e = 35;

% height of plate
m = 256;

% width of plate
n = 256;

% compute number of non-zero values
nz = 2 * m + (n - 2) * (4 * (m - 2) + m);
pz = nz / (n * m)^2;
disp(sprintf("sparsity: %f%%", 100 * pz));

% construct the sparse matrix
t = zeros(nz, 3);
k = 1;
for j = 1 : n
  for i = max((j - 1), 1) : min((j + 1), n)
    for jj = 1 : m
      for ii = max((jj - 1), 1) : min((jj + 1), m)
        if ((i == 1 || i == n || ii == 1 || ii == m) && (i == j) && (ii == jj)) || \
            (((i == j - 1 || i == j + 1)) && (i != 1 && i != n && ii != 1 && ii != m && ii == jj)) || \
            (i == j && (ii == jj - 1 || ii == jj + 1) && ii != 1 && ii != m && i != 1 && i != n)
          t(k, 1) = (i - 1) * m + ii;
          t(k, 2) = (j - 1) * m + jj;
          t(k, 3) = 1;
          k = k + 1;
        else if ((i == j) && (ii == jj) && (ii != 1) && (ii != m))
          t(k, 1) = (i - 1) * m + ii;
          t(k, 2) = (j - 1) * m + jj;
          t(k, 3) = -4;
          k = k + 1;
          endif
        endif
      endfor
    endfor
  endfor
endfor
s = spconvert(t);
disp("constructed matrix");

% solve the eigen system
[v,ls,f] = eigs(s, e, "sm");
l = diag(ls);
disp("solved eigen system");

% save the raw data
fid = fopen(sprintf("modes-%d-%d.dat", m, n), "w", "native");
fwrite(fid, l, "float32", "native");
fwrite(fid, v, "float32", "native");
fclose(fid);
disp("saved raw data");
