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

% sistema
m=100; 
c=50; 
k=2000; 
k3=300; % não linearidade cúbica

% força harmônica
eps=1/3;
F=1500; 
wf=sqrt(k/m)*eps;

up=[u(2); -(c/m)*u(2)-(k/m)*u(1)-(k3/m)*u(1)^3+(F/m)*sin(wf*t)];