Day 4

% Solve a second-order ODE: spring-mass-damper
% m*x'' + c*x' + k*x = F(t)
% State-space form: x1 = x, x2 = x'
% x1' = x2
% x2' = (F(t) - c*x2 - k*x1) / m

m = 1;    % mass (kg)
c = 0.5;  % damping (N·s/m)
k = 4;    % spring constant (N/m)
F = @(t) (t < 1);   % step force: 1 for t<1, then 0

% ODE function: state = [x; x_dot]
ode_fn = @(t, s) [s(2);
                   (F(t) - c*s(2) - k*s(1)) / m];

% Initial conditions: x=0, x'=0
s0   = [0; 0];
tspan = [0, 10];

% Solve with ode45 (Runge-Kutta 4/5 — default for smooth problems)
[t, s] = ode45(ode_fn, tspan, s0);

% For stiff systems, use ode15s instead
% [t, s] = ode15s(ode_fn, tspan, s0);

figure;
subplot(2,1,1);
plot(t, s(:,1), 'b', 'LineWidth', 1.5);
title('Position x(t)'); xlabel('Time (s)'); ylabel('x (m)');
grid on; xline(1, 'r--', 'Force off');

subplot(2,1,2);
plot(t, s(:,2), 'r', 'LineWidth', 1.5);
title('Velocity x''(t)'); xlabel('Time (s)'); ylabel('x'' (m/s)');
grid on;

% Run a Simulink model from MATLAB script
% (assumes 'spring_damper.slx' exists and has Workspace logging)

model = 'spring_damper';
load_system(model);

% Set parameters in the model's workspace
set_param(model, 'StopTime', '10');
set_param(model, 'Solver', 'ode45');
set_param(model, 'MaxStep', '0.01');

% Assign variables used by the model
assignin('base', 'mass', 1.0);
assignin('base', 'damping', 0.5);
assignin('base', 'spring_k', 4.0);

% Run the simulation
simOut = sim(model, 'ReturnWorkspaceOutputs', 'on');

% Access logged signals
t = simOut.get('tout');
x = simOut.get('yout');

figure;
plot(t, x.signals.values);
title('Simulink Output'); xlabel('Time (s)'); ylabel('Position (m)');
grid on;

% Parameter sweep: vary damping coefficient
damping_values = [0.1, 0.5, 1.0, 2.0];
figure; hold on;
for i = 1:length(damping_values)
    assignin('base', 'damping', damping_values(i));
    out = sim(model, 'ReturnWorkspaceOutputs', 'on');
    t_i = out.get('tout');
    x_i = out.get('yout');
    plot(t_i, x_i.signals.values, 'LineWidth', 1.5, ...
         'DisplayName', sprintf('c = %.1f', damping_values(i)));
end
hold off;
legend('Location', 'best');
title('Parameter Sweep: Damping Coefficient');
xlabel('Time (s)'); ylabel('Position (m)');
grid on;

% Create a new Simulink model programmatically
model = 'integrator_demo';
new_system(model);
open_system(model);

% Add blocks
add_block('simulink/Sources/Step',         [model '/Step'],       'Position', [50,50,80,80]);
add_block('simulink/Continuous/Integrator',[model '/Integrator'], 'Position', [150,45,200,85]);
add_block('simulink/Sinks/Scope',          [model '/Scope'],      'Position', [280,47,320,83]);

% Connect blocks: Step -> Integrator -> Scope
add_line(model, 'Step/1',       'Integrator/1');
add_line(model, 'Integrator/1', 'Scope/1');

% Configure step block
set_param([model '/Step'], 'Time', '1', 'After', '1', 'Before', '0');

% Set simulation parameters
set_param(model, 'StopTime', '5', 'Solver', 'ode45');

% Save and run
save_system(model);
sim(model);

% The Scope shows the integral of a step — a ramp