A couple of things going on here. First, it's good to use more points. Use a lot more points. There is no need to go with something as small as 1024 and there is also no need to go with a power of 2, the importance of which, in my opinion, is a myth. I used a nice even million points in the code below.
Second, in the expression for y there is a factor of (eta^2-f.^2) in the denominator. That puts poles at +-eta on the real frequency axis. The path of integration runs right over them, so by implication you are taking the principal value of the integral. The p.v. is the same as the average value of two functions (y3 and y4 below) where the poles are moved slightly off of the real axis in either the plus or minus i direction. The small parameter gamma accomplishes that.
a=0.5;eta=0.3;beta=0.025;wb=14.25;gamma = 1e-10;
plot(f,y),title('Frequency response');ylabel('Response of beam');xlabel('Frequency');
plot(t,(y2)),title('Time response');ylabel('Response of beam');xlabel('Time');
axis([0 2.5 -0.05 0.05]);