Frequency-amplitude and frequency- phase angle solution

The following is the Matlab code for frequency-amplitude and frequency- phase angle solution of mass-spring-damper system with harmonic external excitation force. 

m = 2;                      %  mass of the system
c = 30;                     %  damping of the system
k = 5000;                   %  stiffness of the system

P0 = 1e2;                   %  magnitude of the excitation force
t_final = 3;                %  final time of calculation
n = 5000;                   %  number of data points

wn = sqrt(k/m);             %  natural frequency (rad/sec)
zeta = c/(2*wn*m);          %  viscous damping factor
Input some constants of the system, like mass, stiffness, etc. Calculate the natural frequency of the system and the damping coefficient, zeta. P0, Omega are the amplitude and the excitation frequency of the external harmonic force, repectively. 
Omega = 0;

for i = 1:n
  Omega(i) = 2.5*wn*(i-1)/(n-1);

%
%   Amplitude and phase angle of
%   Particular solution
%

  phi(i) = atan(c*Omega(i)/(k-m*Omega(i)^2));
  D(i) = P0/(sqrt((k-m*Omega(i)^2)^2+(c*Omega(i))^2));

end
The excitation frequency goes from 0 to 2.5 times the natural frequency. We calculated the amplitudes and phase angles for corresponding excitation frequencies. 
iz = find(phi<0);
iz1 = iz(1);

for i = 1:iz1-1
  Phi(i) = phi(i);
end

for i = iz1:n
  Phi(i) = phi(i)+pi;
end
This part is just for nice-looking plot. since 'atan' function in Matlab gives the angel between -pi < angle < pi , you will not get the similar plot that is in the textbook. So, the angles in negative value were added pi (180 degree), which gives the same result. 

%
%   plotting the results
%

figure(1);clf;

subplot(2,1,1),plot(Omega/wn,D/(P0/k));
xlabel('Omega/wn');
ylabel('D/(P0/k)');

subplot(2,1,2),plot(Omega/wn,Phi*180/pi);
xlabel('Omega/wn');
ylabel('phase angle');
axis([0 2.5 0 180]);
Plotting the results. The axis for frequency-phase angle plot was set for the phase angle to show up to 180 degrees.

Click here for the entire Matlab source code for the solution.