function aufg32_6    % Biegeschwingung, links eingespannt, rechts Loslager

EI     = 3000 ;
rhoA   =    3 ;
L      =    1 ;
lamanf =    0 ;   % Nach Nullstellen (Eigenwerten) zu ...
lamend =   80 ;   % ... untersuchender Bereich
nevmax =    3 ;   % Maximale Anzahl zu suchender Eigenwerte

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/L^2) / (2*pi) 

z = 0 : 0.01 : 1 ;
for i = 1:iev-1
    lami = lam0(i) ;
    Z = cos(lami*z)-cosh(lami*z)-(sin(lami*z)-sinh(lami*z))* ...
        (cos(lami)-cosh(lami))/(sin(lami)-sinh(lami)) ;   % Schwingungsformen ...
    subplot (iev-1,1,i) ; plot (z , Z) , grid on , ...    % ... darstellen
             title (strcat(num2str(i),'. Eigenschwingungsform')) ;
end

function f = ewg (lam)
f = sin(lam)*cosh(lam)-cos(lam)*sinh(lam) ;