clc
clear all
close all

%Fenómeno de resonancia

%Ecuación diferencial de segundo orden no homogenea:
% mx"=-k*x-b*x'+F0*sin(w0*t)

%Sistema de ecuaciones diferenciales de primer orden equivalente:
%x'=y
%y'=-(k/m)*x-(b/m)y+(F0/m)*sin(w0*t)

%Parámetros físicos
m=1;% masa
k=1;% constante elástica del resorte
b=0.1;% coeficiente de fricción
F0=1;% fuerza osciladora máxima
x0=1;% posición inicial
y0=0;% velocidad inicial
t0=0;% tiempo inicial
tf=250;% tiempo final

w=(abs((b/(2*m))^2-k/m))^(1/2); % frecuencia angular del sistema
f=w/(2*pi);% frecuencia de la oscilación 
periodo=1/f;% período de la oscilación
w0=0:2*w/40:2*w;% frecuencia de oscilación de la fuerza osciladora

%Método de Runge-Kutta
A=zeros(1,length(w0));
for j=1:length(w0)
h=0.01;% paso
T=t0:h:tf;% vector de tiempos
X=zeros(1,(tf-t0)/h+1);% vector de posiciones
Y=zeros(1,(tf-t0)/h+1);% vector de velocidades
X(1,1)=x0;
Y(1,1)=y0;
for i=2:(tf-t0)/h+1
   k1=Y(1,i-1);
   l1=-(k/m)*X(1,i-1)-(b/m)*Y(1,i-1)+(F0/m)*sin(w0(1,j)*T(1,i-1));
   k2=Y(1,i-1)+h*l1/2;
   l2=-(k/m)*(X(1,i-1)+h*k1/2)-(b/m)*(Y(1,i-1)+h*l1/2)+(F0/m)*sin(w0(1,j)*(T(1,i-1)+h/2));
   k3=Y(1,i-1)+h*l2/2;
   l3=-(k/m)*(X(1,i-1)+h*k2/2)-(b/m)*(Y(1,i-1)+h*l2/2)+(F0/m)*sin(w0(1,j)*(T(1,i-1)+h/2));
   k4=Y(1,i-1)+h*l3;
   l4=-(k/m)*(X(1,i-1)+h*k3)-(b/m)*(Y(1,i-1)+h*l3)+(F0/m)*sin(w0(1,j)*(T(1,i-1)+h));
   X(1,i)=X(1,i-1)+(h/6)*(k1+2*k2+2*k3+k4);
    Y(1,i)=Y(1,i-1)+(h/6)*(l1+2*l2+2*l3+l4);
end
W0=w0./w;
A(1,j)=max(X(1,(100/h):length(X)));
end

figure(1)
plot(W0,A,'+')
title('Amplitud en función de w0/w')
xlabel('w0/w')
ylabel('Amplitud')