% vibration_test
% this is mx'' + cx' + kx = F(t)
% with c = 0, m = 1, k = 4, F(t) = cos(wt)
% resonance occurs at w = 2

w = input('Enter frequency w: ');

options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% [0,20] indicates t range from 0 to 20
% [0,0] is the initial condition x1 = 0, x2 = 0
% Note that ode45 wants a function of t and x as
% first argument so we use an anonymous function
% to make one.

[t,x] = ode45(@(t,x) vibration(t,x,w), [0,20], [0,0], options);

subplot(2,2,1)
plot(t, x(:,1)) % plot just the position x1

subplot(2,2,2)
plot(t, x(:,2)) % plot just the velocity x2

subplot(2,2,3)
plot(t, x) % plot the position and velocity on one plot

subplot(2,2,4)
plot(x(:,1), x(:,2)) % phase plane