clear all
%
b=1;
k=10;
m=10;
%
om0=sqrt(k/m);
gam=b/m;
%
omega=linspace(om0-5*gam,om0+5*gam);
A=om0^2./(om0^2-omega.^2-i*gam*omega);
%
gain =abs(A);
phase=angle(A);
realA=real(A);
imagA=imag(A);
%
omplot=(omega-om0)/gam;
%
subplot(1,2,1)
plot(omplot,gain,'-',omplot,phase,'--','LineWidth',2)
legend('|A|','\phi')
xlabel('(\omega-\omega_0)/\gamma')
hold on
plot([-0.5 -0.5],[0 8],':')
plot([ 0.5  0.5],[0 8],':')
hold off
subplot(1,2,2)
plot(omplot,imagA,'-',omplot,realA,'--','LineWidth',2)
legend('Im(A)','Re(A)')
xlabel('(\omega-\omega_0)/\gamma')
hold on
plot([-5 5],[0 0],':')
hold off
%
print -deps mkb.eps