function two_dof_system()
    % Define our two-dof system
    system.m1 = 0.5;
    system.m2 = 0.5;
    system.k1 = 1;
    system.k2 = 1;
    system.F = 0.1;
    system.time_step = 0.01;
    system.t_final = 50;

    % Explicitly integrate our system
    [t, u1, u2] = integrate_explicit_equations(system);
    
    % Plot the output
    figure(1);
    set(gcf, 'color', 'w');
    clf;
    hold on;
    plot(t, u1, 'DisplayName', 'u1')
    plot(t, u2, 'DisplayName', 'u2');
    box on;
    xlabel('Time');
    ylabel('Displacement');
    legend show;

    % Print the output
    disp(' ');
    disp('*******************')
    disp('SOLUTION:');
    disp(['  average displacement: u1 = ', num2str(mean(u1))]);
    disp(['  average displacement: u2 = ', num2str(mean(u2))]);
    disp('*******************')
end

function [t_out, u1_out, u2_out] = integrate_explicit_equations(system)
    % Unpack some values for convenient access
    dt = system.time_step;
    m1 = system.m1;
    m2 = system.m2;
    k1 = system.k1;
    k2 = system.k2;
    F = system.F;
    
    % Setup initial conditions and time steps
    t_out = 0:dt:system.t_final;
    u1_out = [0];
    u2_out = [0];
    u1dot_out = [0];
    u2dot_out = [0];

    %
    % Explicit time stepping loop
    %
    for t = t_out(2:end)
        % Get end-of-step values for displacement
        u1_n = u1_out(end);
        u2_n = u2_out(end);

        % compute accelerations
        v1dot = (-F - k1*(u1_n - u2_n))/m1;
        v2dot = (-k1*(u2_n-u1_n) - k2*u2_n)/m2;

        % integrate to get velocities
        v1 = u1dot_out(end) + v1dot*dt;
        v2 = u2dot_out(end) + v2dot*dt;

        % integrate to get displacements
        u1 = u1_n + v1*dt;
        u2 = u2_n + v2*dt;

        % save results for output and next step calcs
        u1_out = [u1_out; u1];
        u2_out = [u2_out; u2];
        u1dot_out = [u1dot_out; v1];
        u2dot_out = [u2dot_out; v2];
    end
end

