function aufg32_8det % Biegeschwingung, links eingespannt, Loslager in der Mitte,
                     % Einzelmasse rechts

EI     = 3000 ;
rhoA   =    3 ;
L      =    1 ;
m      =    2 ;

lamanf =    0 ;   % Nach Nullstellen (Eigenwerten) zu ...
lamend =   20 ;   % ... untersuchender Bereich
nevmax =    3 ;   % Maximale Anzahl zu suchender Eigenwerte
LS     =  L/2 ;   % Bezugslaenge L*

global mdm 
mdm   = m/(rhoA*LS) ;

delta = (lamend-lamanf)/1000 ;  % Schrittweite für das "Scannen" des Bereichs
f1    = ewg (lamanf) ;
iev   = 1 ;

for lam = lamanf+delta : delta : lamend
  f2 = ewg (lam) ;
  if (f1*f2 < 0 | f2 == 0)                        % Vorzeichenwechsel? Dann ...
    lam0(iev) = fzero (@ewg , [lam-delta lam]) ;  % ... Nullstelle suchen 
    iev = iev + 1 ;
    if iev > nevmax break ; end                   % nevmax Nullstellen gefunden?
  end
  f1 = f2 ;
end

EigenFrequenzen = (sqrt(3000/rhoA) * lam0.^2/LS^2) / (2*pi) 

zS = 0 : 0.02 : 1 ;
z  = 0 : 0.01 : 1 ;
for i = 1:iev-1
    lami = lam0(i) ;
    c = null (MatrixA(lami)) ; 
    Z( 1: 51) = c(1)*cos(lami*zS)+c(2)*sin(lami*zS)+c(3)*cosh(lami*zS)+c(4)*sinh(lami*zS) ;
    Z(51:101) = c(5)*cos(lami*zS)+c(6)*sin(lami*zS)+c(7)*cosh(lami*zS)+c(8)*sinh(lami*zS) ;
    subplot (iev-1,1,i) ; plot (z , Z) , grid on , ...              % Schwingungsformen ...
             title (strcat(num2str(i),'. Eigenschwingungsform')) ;  % ... darstellen
end

function f = ewg (lam)
f = det(MatrixA(lam)) ;

function A = MatrixA (lam)
global mdm 
A = [  1         0         1         0          0         0         0         0     ;
       0         1         0         1          0         0         0         0     ;
     cos(lam)  sin(lam) cosh(lam) sinh(lam)     0         0         0         0     ;
    -sin(lam)  cos(lam) sinh(lam) cosh(lam)     0        -1         0        -1     ;
    -cos(lam) -sin(lam) cosh(lam) sinh(lam)     1         0        -1         0     ;
       0         0         0         0          1         0         1         0     ;
       0         0         0         0      -cos(lam) -sin(lam) cosh(lam) sinh(lam) ;
       0    0    0    0    sin(lam)+mdm*lam*cos(lam)  -cos(lam)+mdm*lam*sin(lam) ...
                           sinh(lam)+mdm*lam*cosh(lam) cosh(lam)+mdm*lam*sinh(lam)  ] ;