function [] = String_vib_2_L3   % vibration of a string
%......................................................
%      x_1=0  x_2   x_3    x_4    x_5= L
%  T <==*-----(1)----*-----(2)-----*==> T 
% Fixed U_1   U_2   U_3    U_4    U_5 Fixed
% Connectivity :    e     i      j     k
%                   1     1      2     3  % L3
%                   2     3      4     5  % L3
%......................................................
n_e = 2 ; n_g = 1 ; n_n = 3 ; n_m = 5 ;     % constants
n_i = n_n*n_g ; n_d = n_g*n_m ;    % elem & system DOFs
L = 2. ; T = 6e5 ; rho = 0.0234    ; % given m, N, kg/m
X = [0 0.25 0.5 0.75 1]*L          ; % node coordinates
L_e = L / n_e                        ; % element length
Free = [2,3,4]                   ; % free displacements

% Constant quadratic L3 element square matrices
K_e =  T * [ 7, -8,  1, ; ...      %       stiffness L3
            -8, 16, -8, ; ...      %
             1, -8,  7 ] / (3*L_e) ;
M_e = rho *L_e*[4, 2, -1 ; ...     %            mass L3
                2, 16, 2 ; ...     %
               -1,  2, 4 ] / 30    ;

% Assemble two 3x3 element sq matrices into a 5x5
K=zeros (n_d, n_d);  M=zeros (n_d, n_d)    ; % allocate
for k = 1:n_e; % loop over elements ===> ===> ===> ===> 
 rows = [1:n_n] + (k - 1)*(n_n - 1)    ; % connectivity 
 K (rows, rows) = K (rows, rows) + K_e ; % add to stiff
 M (rows, rows) = M (rows, rows) + M_e ;  % add to mass      
end ; % for element k <== <=== <=== <=== <=== <=== <===

%  Solve for eigenvalues and eigenvectors, general form
[Vec, Diag] = eig(K(Free,Free), M(Free,Free))  ; % part
% Vec = eigenvector cols,   Diag = eigenvalue^2 on diag
Omega = sqrt(real(abs(diag(Diag))))  ; % diagonal freqs

fprintf ('String natural frequencies \n')     % heading
exact = pi() * sqrt(T/(rho*L^2))     ; % exact constant
V_abs = max (abs (Vec))          ; % relative amplitude
for k = 1:n_d-2;% loop over DOFs - two EBCs ==> ==> ==>
  omega = Omega(k)                     ; % natural freq
  true  = k * exact              ; % exact natural freq
  fprintf ('%i, FEA = %8.3e, Exact = % 8.3e \n', ...
            k, omega, true)              % pretty print
  Vec (:, k) = Vec (:, k) / V_abs (k) ; % scale abs = 1
end  ; % for all active DOF <== <== <== <== <== <== <==

%    Pad mode shape with end zero BC displacements
Pad = [0, 0, 0; Vec; 0, 0, 0]   ; % insert rows of zero
fprintf ('          Mode Shapes \n')        ; % heading
fprintf ('Node     1         2         3 \n');% heading
for k = 1:n_m   ; % loop over nodes ---> ---> ---> --->
  fprintf ('%i    %9.3e  %9.3e  %9.3e \n', ...
            k, Pad(k,:))             ; % relative shape
end   ; % for node k relative shape <--- <--- <--- <---
% end of String_vib_2_L3 ==============================
% Running gives:
% String natural frequencies 
% 1, FEA = 7.984e+03, Exact =  7.954e+03 
% 2, FEA = 1.601e+04, Exact =  1.591e+04 
% 3, FEA = 2.873e+04, Exact =  2.386e+04
%          Mode Shapes
% Node     1          2          3
% 1    0.000e+00  0.000e+00  0.000e+00
% 2    7.068e-01 -1.000e+00  4.068e-01
% 3    1.000e+00  1.473e-16 -1.000e+00
% 4    7.068e-01  1.000e+00  4.068e-01
% 5    0.000e+00  0.000e+00  0.000e+00