How to obtain the frequency response of a linear Multi-Degree of Freedom system

Hi all,

I have the following Multi-Degree of Freedom (MDOF) system for which I would like to find its frequency response:

where and U is the step input and is the first undamped natural frequency.

Being a simple second order system it seems pretty straight forward to find its frequency response to the unit load, however when i write the code in MATLAB, I get a very different result from the solution of the excercise (see below). This is even more puzzling as the solution is available here (pg. 26) and I believe I am performing the same process.

So I would like to ask, is there something that I am forgetting or getting wrong within MATLAB from what you can see? I think something wrong might be happening in the inversion of the (s^2*eye(n) + s*(V'*C*V) + D) matrix (maybe numerical problems?)

clear ; clc

% Input
n = 3 ;         % Number of DOFs
m = 1 ;         % Mass [kg]
k = 1 ;         % Stiffnes [N/m]
c = 0.2 ;       % Viscous damping [Ns/m]
nv = [1 0 0]' ; % Load application vector (applied to 1st DOF)

% Matrices
M = diag(m*ones(n, 1)) ;
K = diag(2*k*ones(n,1)) - diag(k*ones(n-1,1),-1) - diag(k*ones(n-1,1),1) ;
C = diag(c*ones(n,1)) ;

% Undamped Natural Frequencies & Modes
[V, D] = eig(K, M) ;
Om = diag(sqrt(abs(D))) ; % Natural frequencies vector

% Frequency Aanalysis
t0 = 2*pi/Om(1) ;           % Time where unit step is applied
omega = linspace(1e-3, 3) ; % Frequency vector

for i = 1:numel(omega)
    
    s = 1i*omega(i) ;
    
    Us = V'*((1 - exp(-s*t0)/s)*nv) ;
    FRF(:,i) = (s^2*eye(n) + s*(V'*C*V) + D)\Us ;
end

% Plotting
figure ; plot(omega, abs(FRF)) ; ylim([0 10]) ; grid on ; set(gca, 'YScale', 'log') ;

[Danillo Gaspar]

Have you found a solution?