Can someone help me to obtain the FRF curve and plot the curve?

%%
% clear memory
clear all
clc

%% % Define material properties and beam dimensions % % %
%xx—————————————————————————————-xx%
width = 0.3; % beam width
t1 = 0.1; % beam thickness of first layer
t2 = 0.1; % beam thickness of second layer
L = 1; % beam length
poisson1 = 0.33; % Poisson’s ratio first layer
poisson2 = 0.33; % Poisson’s ratio second layer
rho1 = 1700; % density of material 1
rho2 = 1240; % density of material 2
E1 = 2.2e11; % modulus of elasticity for material 1
E2 = 7.9e19; % modulus of elasticity for material 1
E3=7.9e9; % modulus of elasticity for material 2

% % % beam load parameters % % %
%xx—————————————————————————————-xx%
n = 0; % material gradation index
P = -100; % uniform distributed load
Mb = 24.4; % ratio of Young’s moduli between two materials
E2 = E1 / Mb;

% % % Derived FGM parameters % % %
%xx—————————————————————————————-xx%
[C1, I01, I1, A01] = effectstifmass(width, t1, rho1, rho2, E1, E2, Mb, n, poisson1);
[C2, I02, I2, A02] = effectstifmass(width, t2, rho1, rho2, E3, E3, Mb, n, poisson2);

% Mass parameters for FGM
numberElements = 40;
nodeCoordinates = linspace(0, L, numberElements + 1);
xx = nodeCoordinates’;
x = xx’;
elementNodes = zeros(size(nodeCoordinates, 2) – 1, 2);
for i = 1:size(nodeCoordinates, 2) – 1
elementNodes(i, 1) = i;
elementNodes(i, 2) = i + 1;
end
numberNodes = size(xx, 1);
GDof = 2 * numberNodes;

% % % Static analysis % % %
%xx—————————————————————————————-xx%
% Computation of the system stiffness, force, mass
[stiffness1, force, mass1] = formStiffnessMassTimoshenkoBeam(GDof, numberElements, …
elementNodes, numberNodes, xx, C1, P, rho1, I01, t1);
[stiffness2, ~, mass2] = formStiffnessMassTimoshenkoBeam(GDof, numberElements, …
elementNodes, numberNodes, xx, C2, P, rho1, I02, t2);

% % % Boundary conditions % % %
%xx—————————————————————————————-xx%
fixedNodeW = [1, 1];
fixedNodeTX = [1, 2];
prescribedDof = [fixedNodeW; fixedNodeTX + numberNodes];
stiffness = stiffness1 + stiffness2 + stiffness1;
mass = mass1 + mass2 + mass1;
K = stiffness;
m = mass;

% % % Free vibration analysis % % %
%xx—————————————————————————————-xx%
[modes, eigenvalues] = eigenvalue(GDof, prescribedDof, stiffness, mass, 0);
omega = sqrt(eigenvalues) * L * L * sqrt(rho1 * A01 / E1 / I01);
disp(‘First 3 dimensionless frequencies:’)
disp(omega(1:3))
f = sqrt(eigenvalues);
disp(‘First 3 natural frequencies (Hz):’)
disp(f(1:3))

% Drawing mesh and deformed shape
modeNumber = 3;
V1 = modes(:, 1:modeNumber);

% Normalize the mode shapes for better visualization
for i = 1:modeNumber
V1(:,i) = V1(:,i) / max(abs(V1(:,i)));
end

% Plot the first 3 mode shapes
figure;
for i = 1:modeNumber
subplot(3,1,i);
plot(x, V1(1:numberNodes, i), ‘LineWidth’, 2);
title([‘Mode Shape ‘, num2str(i)]);
xlabel(‘Beam Length (x)’);
ylabel(‘Normalized Displacement’);
grid on;
end

% Solve for static displacements
displacements = solution(GDof, prescribedDof, stiffness, force);
U = displacements;

% Max displacement
disp(‘Max normal displacement:’)
disp(min(U(1:numberNodes)))

% % % Forced vibration analysis % % %
%xx—————————————————————————————-xx%
c1 = 1.5;
c2 = 0.9; % damping coefficients
Ca = c1 * mass + c2 * stiffness;
F = @(t) force * sin(2* t); % time-dependent load
D0 = zeros(GDof,1); % initial displacement
V0 = zeros(GDof,1); % initial velocity
dt = 0.005; T = 100; tt = 0:dt:T;
[D, V, A] = Newmark(mass, Ca, stiffness, F, D0, V0, dt, T);

% Plot forced vibration displacement at a node
figure;
plot(tt, D(21,:), ‘b’, ‘LineWidth’, 2);
title(‘Forced Vibration Displacement at Node 21’);
xlabel(‘Time (s)’);
ylabel(‘Displacement’);
grid on;

% % % Export results to Excel % % %
%xx—————————————————————————————-xx%
xlswrite(‘mode vibration result.xlsx’, [x(:), V1(1:numberNodes,:)]);
xlswrite(‘normal displacement.xlsx’, [x(:), U(1:numberNodes,:)]);
xlswrite(‘forced displacement.xlsx’, [tt’, D(21,:)’]);
this is my current code%%
% clear memory
clear all
clc