Calculating FRF for MDoF Euler-Bernoulli beam using specified formula

I'm dealing with a Euler-bernoli beam which has 40 elements (ne = 40) and 80 Degrees of Freedom (DoFs) (nd = 80). I've calculated Mode Shapes (V) and Natural frequencies (D) using eig built-in function. Now I need to find the FRF (Frequency Response Function) for my data using the FRF formula for MDoF systems:

in which k is the modeshape number (from 1 to N which is 80 in my case) and is the i-th degree of freedom of k-th modeshape, likewise is the j-th degree of freedom of k-th modeshape (correct me if I'm wrong)
is the natural frequeny of the k-th modeshape:
[w,index]=sort(abs(diag(D)));
is the frequency of the applied load (which I have 80 and I'm planning on using each frequency on the corresponding modeshape) I'm using w2 to represent this operand like below:
w2=(w.^0.5)/(2*pi);
My structure's DoFs for 40 elements are like below which first row denotes start dofs and 2nd row denotes end dofs for all the 40 elements:
dof=[1:2:79; 4:2:80,80]
I need to find which should be 80x80