function Aufg32_9     % Biegeschwinger (Aufgabe 32-9)

EI   = 5000 ;
rhoA =  4.8 ;
cT   =  60000 ;
a    =  0.6 ;
b    =  0.8 ;
m    =    4 ;

lamanf =    0 ;   % Nach Nullstellen (Eigenwerten) zu ...
lamend =   20 ;   % ... untersuchender Bereich
nevmax =    3 ;   % Maximale Anzahl zu suchender Eigenwerte

global mdrAa bda cTadEI
mdrAa  = m/(rhoA*a) ;
bda    = b/a        ;
cTadEI = cT*a/EI   ;

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(EI/rhoA) * lam0.^2/a^2) / (2*pi) 

dz   = (a+b)/1400      ;
z    =   0 : dz : a+b  ;
z1da = ( 0 : dz : a)/a ;
z2da = (dz : dz : b)/a ;

for i = 1:iev-1
    lami = lam0(i) ;
    A = [   1          0          1          0            0            0              0        ;
           lami      cTadEI     -lami     cTadEI          0            0              0        ;
         cos(lami)  sin(lami) cosh(lami) sinh(lami)      -1            0             -1        ;
        -sin(lami)  cos(lami) sinh(lami) cosh(lami)       0           -1              0        ;
        -cos(lami) -sin(lami) cosh(lami) sinh(lami)       1            0             -1        ;
         sin(lami) -cos(lami) sinh(lami) cosh(lami)   mdrAa*lami       1          mdrAa*lami   ;
            0          0          0          0      cos(bda*lami) sin(bda*lami) cosh(bda*lami) ] ;
    B = [    0     ;
             0     ;
             0     ;
             1     ;
             0     ;
             1     ;
        -sinh(bda*lami)] ;
    c = inv(A)*B ;
    Z1 = c(1)*cos(lami*z1da)+c(2)*sin(lami*z1da)+c(3)*cosh(lami*z1da)+c(4)*sinh(lami*z1da) ;
    Z2 = c(5)*cos(lami*z2da)+c(6)*sin(lami*z2da)+c(7)*cosh(lami*z2da)+     sinh(lami*z2da) ;
    subplot (iev-1,1,i) ; plot (z , [Z1 Z2]) , grid on , ...        % Schwingungsformen ...
             title (strcat(num2str(i),'. Eigenschwingungsform')) ;  % ... darstellen
end

function f = ewg (lam)
global mdrAa bda cTadEI
A = [   1         0         1         0           0           0              0            0        ;
       lam     cTadEI     -lam     cTadEI         0           0              0            0        ;
     cos(lam)  sin(lam) cosh(lam) sinh(lam)      -1           0             -1            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        ;
     sin(lam) -cos(lam) sinh(lam) cosh(lam)   mdrAa*lam       1          mdrAa*lam       -1        ;
        0         0         0         0      cos(bda*lam) sin(bda*lam) cosh(bda*lam) sinh(bda*lam) ;
        0         0         0         0     -sin(bda*lam) cos(bda*lam) sinh(bda*lam) cosh(bda*lam) ] ;
f = det(A) ;