Try this code:
clear; clc;
% Define beam parameters
L = 2; % Length of the beam
num_elements = 100; % Number of elements
num_modes = 6; % Number of modes to calculate
% Prompt user for beam properties
rho = 2711.518;
A = 0.00555;
I = 0.00003186625;
E = 71 * 10^9;
% Compute element size
dx = L / num_elements;
% Assemble mass and stiffness matrices
M = zeros((num_elements+1) * 2, (num_elements+1) * 2);
K = zeros((num_elements+1) * 2, (num_elements+1) * 2);
for i = 1:num_elements
% Element mass matrix
m = [156 22*dx 54 -13*dx;
22*dx 4*dx^2 13*dx -3*dx^2;
54 13*dx 156 -22*dx;
-13*dx -3*dx^2 -22*dx 4*dx^2] * (rho * A * dx / 420);
% Element stiffness matrix
k = [12 6*dx -12 6*dx;
6*dx 4*dx^2 -6*dx 2*dx^2;
-12 -6*dx 12 -6*dx;
6*dx 2*dx^2 -6*dx 4*dx^2] * (E * I / dx^3);
% Assemble into global mass and stiffness matrices
idx = (i-1)*2 + 1;
M(idx:idx+3, idx:idx+3) = M(idx:idx+3, idx:idx+3) + m;
K(idx:idx+3, idx:idx+3) = K(idx:idx+3, idx:idx+3) + k;
end
% Apply boundary conditions (cantilever beam)
M = M(3:end, 3:end);
K = K(3:end, 3:end);
% Solve eigenvalue problem
[V, D] = eig(K, M);
% Sort eigenvalues and eigenvectors
[D, sort_idx] = sort(diag(D));
V = V(:, sort_idx);
% Extract first num_modes eigenvalues and eigenvectors
omega = sqrt(D(1:num_modes));
modes = V(:, 1:num_modes);
% Plot x-y axis
x = linspace(0, L, (num_elements+1) * 2);
y = zeros(size(x));
figure;
hold on;
title('Cantilever Beam Mode Shapes');
xlabel('Beam Length');
ylabel('Displacement');
plot(x, y, 'k-');
axis equal;
grid on;
% Display natural frequencies
disp('Natural Frequencies:');
disp(omega);
