% Differenzenverfahren, Biegeschwingungen (Beispiel 3)
clear all
L    = 1    ;
rhoA = 3    ;
EI   = 3000 ;
m    = 2    ;

nA  = 4000 ;            % Anzahl der Abschnitte
h   = L/nA ;
n   = nA   ;
A   = zeros (n,3) ;     % Nullmatrix
B   = ones  (n,1) ;     % Einheitsmatrix

A(1:2,:) = [7 -4  1 ; 6 -4  1] ;
for i=3:n-2
    A(i,:) = [6 -4  1] ;
end
A(n-1:n,:) = [5 -2 0 ; 1 0 0] ;

nMitte        = round(n/2+0.1) ;
A(nMitte-2,3) = 0 ;
A(nMitte-1,2) = 0 ;
A(nMitte  ,:) = [1 0 0] ;
B(nMitte  ,1) = 0 ;
B(n       ,1) = 0.5 + m/(rhoA*h) ;

[nc ev Z] = isiasb_m (A,B,3) ;
Eigenfrequenzen = sqrt(ev(1:nc)*EI/(rhoA*h^4))/(2*pi)

z = 0 : h : L ;
for i = 1:nc
    subplot (nc,1,i) ; plot (z , [0 ; Z(:,i)]) , grid on , ... 
             title (strcat(num2str(i),'. Eigenschwingungsform')) ;
end