function dvdt = flinsys(t,v,p) x = v(1); y = v(2); omega = p(1); a = p(2); b = p(3); k = p(4); dvdt = [ y; -b*y - k*x + a*cos(omega*t)];