clear; clc; close all;

% --- 1. System Parameters ---
reg_no = 6562;

p.m1 = reg_no/10;   % Mass 1 (kg)
p.m2 = reg_no/20;   % Mass 2 (kg)
p.c1 = reg_no/10;   % Damping 1 (Ns/m)
p.c2 = reg_no/20;   % Damping 2 (Ns/m)
p.k1 = reg_no;      % Stiffness 1 (N/m)
p.k2 = reg_no/2;    % Stiffness 2 (N/m)
p.F  = (reg_no/100) + 0.50;  % Force Amplitude (N)
p.f  = reg_no/1000;     % Force Frequency (Hz)

% --- 2. Simulation Settings ---
tspan = [0 20];         % Time range (s)
y0 = [0; 0; 0; 0];      % Initial Conditions: [x1; x2; v1; v2]

% --- 3. Define Dynamics System directly ---
% y(1)=x1, y(2)=x2, y(3)=v1, y(4)=v2
sys_dynamics = @(t, y) [
    y(3); ... 
    y(4); ... 
    (-(p.k1 + p.k2)*y(1) + p.k2*y(2) - (p.c1 + p.c2)*y(3) + p.c2*y(4) + p.F * sin(2 * pi * p.f * t)) / p.m1; ... 
    (p.k2*y(1) - p.k2*y(2) + p.c2*y(3) - p.c2*y(4)) / p.m2 ... 
];

% --- 4. Run Solver ---
[t, y] = ode45(sys_dynamics, tspan, y0);

% --- 5. Visualization ---
% Octave on mycompiler supports standard subplots cleanly!
subplot(2,1,1);
plot(t, y(:,1), 'b-', 'LineWidth', 2); hold on;
plot(t, y(:,2), 'r-', 'LineWidth', 2);
ylabel('Displacement (m)');
xlabel('Time (s)');
title('Displacement Response');
legend('Mass 1 (x1)', 'Mass 2 (x2)');
grid on;

subplot(2,1,2);
plot(t, y(:,3), 'b-', 'LineWidth', 2); hold on;
plot(t, y(:,4), 'r-', 'LineWidth', 2);
ylabel('Velocity (m/s)');
xlabel('Time (s)');
title('Velocity Response');
legend('Mass 1 (v1)', 'Mass 2 (v2)');
grid on;