function up=leh2(t,u)
% sistema linear, excitação harmônica, c.i. não nulas, 2 gdl

% sistema
M=[2 0; 0 1]; 
C=[3 -0.5; -0.5 0.5]; 
K=[3 -1; -1 1];
A=[zeros(2) eye(2); -inv(M)*K -inv(M)*C];

% excitação
F=[1; 1]; 
wf=2;
f=inv(M)*F;

up=A*u+[0; 0; f]*sin(wf*t);