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 + cos(2*t);

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


% 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");