function du = duffDE(t,u)
% Duffing's ODE
ep = 0.2;
du1 = u(2);
du2 = -u(1) - ep*(u(1))^3;
du = [du1;du2];
end

