%resonance02.m
%Plot Amplitude

w1=[0:0.001:sqrt(1/2)-0.001];
w2=[sqrt(1/2)+0.001:0.001:sqrt(3)-0.001];
w3=[sqrt(3)+0.001:0.001:4];
w=[w1 w2 w3];

fo=1;
k=1;
m=1;

a1=(3*fo/5)*(1./(5/9*k-10/9*m*w1.^2)+1./(5*k-5/3*m*w1.^2));
a2=(3*fo/5)*(1./(5/9*k-10/9*m*w2.^2)+1./(5*k-5/3*m*w2.^2));
a3=(3*fo/5)*(1./(5/9*k-10/9*m*w3.^2)+1./(5*k-5/3*m*w3.^2));

b1=(-3*fo/5)*(1./(5/9*k-10/9*m*w1.^2)-(2/3)./(5*k-5/3*m*w1.^2));
b2=(-3*fo/5)*(1./(5/9*k-10/9*m*w2.^2)-(2/3)./(5*k-5/3*m*w2.^2));
b3=(-3*fo/5)*(1./(5/9*k-10/9*m*w3.^2)-(2/3)./(5*k-5/3*m*w3.^2));

figure(1)
plot(w1,a1,'r-',w2,a2,'r-',w3,a3,'r-')
axis([0 2 -20 20])
xlabel('\omega')
ylabel('a_{1}')
title('Amplitude vs. Frequency')
hold on
plot([1/sqrt(2) 1/sqrt(2)],[-20,20],'k:')
plot([sqrt(3) sqrt(3)],[-20,20],'k:')
plot([0 4],[0,0],'k-')
hold off

figure(2)
plot(w1,b1,'r-',w2,b2,'r-',w3,b3,'r-')
axis([0 2 -20 20])
xlabel('\omega')
ylabel('a_{2}')
title('Amplitude vs. Frequency')
hold on
plot([1/sqrt(2) 1/sqrt(2)],[-20,20],'k:')
plot([sqrt(3) sqrt(3)],[-20,20],'k:')
plot([0 4],[0,0],'k-')
hold off
