function xprime = vibration(t,x,w)
% this is mx'' + cx' + kx = F(t)
% with c = 0, m = 1, k = 4, F(t) = cos(wt)
% Use @(t,x) vibration(t,x,w) in ode45

xprime(1) = x(2);
xprime(2) = cos(w*t) - 4*x(1);
xprime= [xprime(1); xprime(2)]; % column vector
end