function ex16()
m = 1;
c = 1;

M = m*eye(8,8);         % Massenmatrix
C = c* ...              % Steifigkeitsmatrix
   [2.5,  0.5, -1.0,    0, -0.5, -0.5,    0,    0; ...
    0.5,  1.5,    0,    0, -0.5, -0.5,    0, -1.0; ...
   -1.0,    0,  2.0,    0,    0,    0,    0,    0; ...
      0,    0,    0,  1.0,    0, -1.0,    0,    0; ...
   -0.5, -0.5,    0,    0,  2.0,    0, -1.0,    0; ...
   -0.5, -0.5,    0, -1.0,    0,  2.0,    0,    0; ...
      0,    0,    0,    0, -1.0,    0,  1.5,  0.5; ...
      0, -1.0,    0,    0,    0,    0,  0.5,  1.5];

[U, D] = eig(C,M);
f = sqrt(diag(D))/(2*pi);

fprintf("Eigenfrequenzen")
f
fprintf("Eigenschwingungen")
U

plotMode(U(:,1));
F = getframe(gcf);
imwrite(F.cdata, "bild51.png");

plotMode(U(:,2));
F = getframe(gcf);
imwrite(F.cdata, "bild52.png");

%------------------------------------------------------------
function plotMode(xe)
% plottet das Modell sowie die Eigenschwingung xe

x0 = [1 2 2 1 0 3; 0 0 1 1 0 0];  % Koordinaten der Punkte im Gleichgewicht
A = [0 1 1 1 1 0; 1 0 1 0 0 1; 1 1 0 1 0 1; ...
     1 0 1 0 1 0; 1 0 0 1 0 0; 0 1 1 0 0 0]; % Verbindungen

xp = x0(1,:)';
yp = x0(2,:)';
M = length(xp);
N = length(xe)/2;

% plot masses
plot(xp, yp, "ko")
hold("on")

% plot springs
for I=1:M
  for J = 1:(I-1)
    if A(I,J) ~= 0
       plot([xp(I) xp(J)], [yp(I) yp(J)], "k-")
    end
  end
end
   
% plot vectors
for I=1:N
  quiver(xp(1:N), yp(1:N), xe(1:2:end), xe(2:2:end), "r-", ...
         "AutoScaleFactor",0.5);
end
axis("equal")
axis([-0.5 3.5 -0.5 1.5]);
axis("off")
hold("off")
