In the present paper an elementary symbolic method for the estimation of basins of attraction in second order nonlinear dynamical systems is formulated and its implementation using Mathematica® is shown. The estimation algorithm is based upon the construction of positively invariant compact boxes that trap the trajectories with unbounded initial conditions. We obtain such boxes through a Lyapunov function whose orbital derivative is bounded by a bounding function that can be represented as the addition of two scalar functions. The detection of persistent oscillating behaviors deserves a prominent place between all possible applications of our tool, as it will be shown by considering Fitzhugh Equations as a case study.