TWO DOF SYSTEM EQUATION: PLOTTING FRF(realpart vs frequency ratio)

wn = input('natural frequency of main system, wn:'); z = input('damping factor,z:'); ms = input('mass of main system, ms:'); md = input('mass of damper system, md:'); % Tuned damper parameters kd = input('Stiffness of damper, kd:'); zd = input('damping factor, zd:'); wd = sqrt(kd/md); w = 0:1:1000; % Non-dimensional parameters r = w/wn; rd = wd/wn; m = md/ms; % Fixed for a setup for s=1: length(r) g(s,:) = ((1-r.^2-r.^2).* m.*(((rd.^2./r.^2)-1+ 4*zd.^2)./(((rd./r)-(r./rd)).^2 +4*zd.^2))); h(s,:) = (2*r*z + r.^2).* m.*(((2*zd*r./rd)./(((rd./r)-(r./rd)).^2 +4*zd.^2))); G0(s,:)= g(s,:)./((g(s,:)).^2+(h(s,:)).^2); end figure (1) plot(r,G0) xlabel ('\omega/\omega_n') ylabel ('|x (i\omega)|/F')

I have written this code for plotting FRF of 2DOF system eq. The plot that i am getting is not matching with original FRF. Kindly help me in getting the graph. /correct any mistakes in coding