% Use ode45 from t = 0 to t = 10. Initial conditions are
% x(0) = 1, x'(0) = 0
[t,ysol] = ode45(@rhs,[0,10],[1;0]);

% Make plot in one go
plot(t,ysol); % Plots both columns of 'ysol' in one go
hold on;
plot(t,cos(2*t));


% Make plot in n steps, where n = number of state variables
% plot(t,ysol(:,1)); % Plot first column of ysol
% hold on;
% plot(t,ysol(:,2)); % Plot second column of ysol
% plot(t,cos(2*t));

legend("x","x dot","forcing","Location","southeast");

function dydt = rhs(t,y)
% Implements the right-hand-side of the equation
% dy/dt = f(y,t)
% where y is a vector containing n state variables.
% Here, we are interested in the equation m x'' + k x = cos(2t)
% So y1 = x and y2 = x'.
    % un-pack contents of y
    y1 = y(1);
    y2 = y(2);

    % define constants
    k = 5; m = 3;

    % Define what dy/dt is
    dy1dt = y2;
    dy2dt = -(k/m)*y1 + 5*f2(t);

    % Assemble contents of dy/dt into a column vector
    dydt = [dy1dt;dy2dt];
end

function output = f2(t)
    % Builds a 'hat' centered at t = 2.5
    % using ramp functions.
    output = + heaviside(t-2).*5.*(t-2) + ...
             - heaviside(t-2.2).*5.*(t-2.2) + ...
             - heaviside(t-2.8).*5.*(t-2.8) + ...
             + heaviside(t-3).*5.*(t-3);
end