function xdot = twomassrhs(t,x)
k = 1;
m = 1;
omega = sqrt(3*k/m)-0.08;
F = 10;
A = [0 1 0 0;
    -2*k/m 0 k/m 0;
    0 0 0 1;
    k/m 0 -2*k/m 0];

xdot = A*x;
xdot(4) = xdot(4) + F/m*cos(omega*t);
end