clc
clear all
close all

%Oscilación forzada

%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;% Módulo máximo de la fuerza osciladora
w0=1;% Frecuencia de oscilación de la fuerza osciladora
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

%Método de Runge-Kutta
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:length(T)
   k1=Y(1,i-1);
   l1=-(k/m)*X(1,i-1)-(b/m)*Y(1,i-1)+(F0/m)*sin(w0*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*(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*(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*(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

figure(1)
plot(T,X,'b')
title('Posición en función del tiempo')
xlabel('Tiempo')
ylabel('Posición')

figure(2)
plot(T,Y,'b')
title('Velocidad en función del tiempo')
xlabel('Tiempo')
ylabel('Velocidad')

figure(3)
plot(X,Y,'r',X(1,1),Y(1,1),'b*',X(1,(tf-t0)/h+1),Y(1,(tf-t0)/h+1),'bo')
title('Diagrama de fases')
xlabel('Posición')
ylabel('Velocidad')
legend('X vs V','Estado inicial','Estado final')