clc;
clear;

% 1. ECUACION DIFERENCIAL DE MOVIMIENTO
% 0.5*x'' + 4*x' + 32*x = 10*cos(4*t)
% x'' + 8*x' + 64*x = 20*cos(4*t)

% 2. CONVERSION A SISTEMA DE PRIMER ORDEN
% x1 = x
% x2 = x'
% x1' = x2
% x2' = 20*cos(4*t) - 8*x2 - 64*x1

function dx=sistema(t,x)
    dx(1) = x(2);
    dx(2) = 20*cos(4*t)-8*x(2)-64*x(1);
    dx = [dx(1),dx(2)];
endfunction

% 3. Solucion Numerica
I = [0:0.05:5];

% Condiciones iniciales
x0 = [0.5,-2];

% INGRESAR ODE45
[t,x] = ode45(@sistema,I,x0);

% 4. Resultados
disp("Tiempo, Posicion y Velocidad");
disp([t,x]);

% 6. Analisis

% Primer paso por la posicion de equilibrio
i = find(x(1:end-1,1).*x(2:end,1)<=0,1);

t_equilibrio = t(i)-x(i,1)*(t(i+1)-t(i))/(x(i+1,1)-x(i,1));

% Amplitud para tiempos mayores o iguales a 3 segundos
indices = find(t>=3);

amplitud = (max(x(indices,1))-min(x(indices,1)))/2;

disp("Primer paso por el equilibrio:");
disp(t_equilibrio);

disp("Amplitud aproximada:");
disp(amplitud);

% 5. GRAFICA
plot(t,x(:,1),"-o");

title("Movimiento de la masa amortiguada");
xlabel("Tiempo (s)");
ylabel("Posicion (ft)");
grid("on");
legend("x(t)");