% Aufgabe 32-9 mit Differenzenverfahren
clear all
EI   =   5000 ;
rhoA =    4.8 ;
cT   =  60000 ;
a    =    0.6 ;
b    =    0.8 ;
m    =      4 ;
L    =    a+b ;

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

A(1:2,:) = [6-(1-cT*h/(2*EI))/(1+cT*h/(2*EI)) -4  1 ; 6 -4  1] ;
for i=3:n-2
    A(i,:) = [6 -4  1] ;
end
A(n-1:n,:) = [6 -4 0 ; 7 0 0] ;

nMasse = round (nA*a/L + 0.1) ;
B(nMasse,1) = 1 + m/(rhoA*h)  ;      % Einzelmasse m

[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) ; 0]) , grid on , ... 
             title (strcat(num2str(i),'. Eigenschwingungsform')) ;
end