clc
clear all
close all

w0=1;
F0=1;
m=1;
vector_b=[0.2 0.5 0.8 1];
w=[0:0.01:2*w0];
matriz_A=zeros(length(w),length(vector_b));
vector_wAmax=zeros(1,length(vector_b));
vector_Amax=zeros(1,length(vector_b));
for i=1:length(vector_b)
    b=vector_b(i);
matriz_A(:,i)=(F0/m)*((w.^2-w0).^2+(b*w/m).^2).^(-1/2);
vector_wAmax(1,i)=(w0^2-(1/2)*(b/m)^2)^(1/2);
vector_Amax(1,i)=F0/(b*(w0^2-(b/(2*m))^2)^(1/2));
end
figure(1)
plot(w,matriz_A,'-b',vector_wAmax,vector_Amax,'*k')
hold on
ylabel('Amplitud (m)')
xlabel('\omega (rad/s)')
% La curva que pasa por todos los máximos
d=find(w<w0);
cA=(F0/m)*(w0^4-w(1:max(d)).^4).^(-1/2);
plot(w(1:max(d)),cA,'--k')