If you use a for loop, you need to index X and phi otherwise you will of course overwrite the previous result. So you could do something like:
r = 0:0.1:2.5;
X = zeros(size(all_r));
phi = zeros(size(all_r));
for idx = 1:numel(r)
X(idx) = (f_0/omega_n^2)/(sqrt((1-r(idx)^2)^2+(2*zeta*r(idx))^2)); %Amplification Equation
phi(idx) = atan2((2*zeta*r(idx)),(1-r(idx)^2)); %Phase angle equation
end
But of course, matlab can operate on whole vectors at once so a loop is not needed at all:
r = 0:0.1:2.5;
X = (f_0/omega_n^2)./(sqrt((1-r.^2).^2+(2*zeta*r).^2)); %Amplification Equation
phi = atan2((2*zeta*r),(1-r.^2)); %Phase angle equation