l = 20;

dom = zeros(l, l);
dom(2:l-1, 2:l-1) = ones(l-2, l-2);
dom(l/2+1:l, 1:l/2) = zeros(l/2, l/2);

idx = dom;

pnt = zeros(sum(sum(dom)),2);

n = 0;
for j = 1:l
  for k = 1:l
    if(dom(j,k)>0)
      n++;
      idx(j,k) = n;
      pnt(n,:) = [j,k];
    endif
  endfor
endfor

lap = -4*eye(n,n);

nb = [0,1; 0,-1; 1,0; -1,0];
for j=1:n
  for k=1:4
    nn = pnt(j,:) + nb(k,:);
    if(dom(nn(1),nn(2)) > 0)
      lap(j,idx(nn(1),nn(2))) = 1;
    endif
  endfor
endfor

[ev,ew] = eig(lap);
omega = sqrt(-diag(ew));

elo0 = zeros(l, l);
[x,y] = meshgrid(linspace(-pi/2,pi/2,l/2),linspace(-pi/2,pi/2,l/2));
# elo0(1:l/2, l/2+1:l) = 20*cos(x).*cos(y);
elo0(1:l/2, 1:l/2) = 20*cos(x).*cos(y);

u0 = zeros(n,1);
for j=1:n
  u0(j) = elo0(pnt(j,1),pnt(j,2));
endfor

coeff = ev'*u0;

for t=0:0.5:80
  ut = ev*(cos(omega*t).*coeff);
  for j=1:n
    elo(pnt(j,1),pnt(j,2)) = ut(j);
  endfor
  surf(elo);
  axis([0,20,0,20,-20,20]);
  replot
endfor

k = 0;
for t=0:0.5:80
  ut = ev*(cos(omega*t).*coeff);
  for j=1:n
    elo(pnt(j,1),pnt(j,2)) = ut(j);
  endfor
  vtk_clear;
  vtk_axis([0,20,0,20,-20,20]);
  vtk_surf(elo);
  vtk_axis("off");
  name = sprintf("bild%03i.jpg",k);
  vtk_print(name,"-djpeg");
  k++;
endfor
