function up=aqeh(t,u)
% sistema não linear, amort. quadrático, excitação harmônica (c.i. nulas)

% sistema
m=10; 
af=0.005; % amort. quadrático
k=200; 

% força harmônica
Fo=150; 
wf=sqrt(k/m);

up=[u(2); -(af/m)*sign(u(2))*u(2)^2-(k/m)*u(1)+(Fo/m)*sin(wf*t)];