function s = MDOF_simulation(t,y,flag,M,C,K)
%
s = zeros(4,1);
A=[zeros(2,2) eye(2,2);-M\K -M\C];
s=A*y;