Chaotic motions of a two dimensional airfoil with coupled structural nonlinearities, both in pitch as well as plunge degrees of freedom, are investigated via a numerical integration method. The original system of coupled integro-differential equations governing the motion of the present aeroelastic model is transformed into a simple system of six ordinary differential equations (ODEs), rather than the previously frequently used eight ODEs. Complex dynamical behaviors are revealed and identified through the means of bifurcation diagrams, the phase portraits, the amplitude spectra and the Poincare maps. Besides, a more quantitative method, namely that of observing the evolution of the largest Lyapunov exponent (LLE) is also applied to diagnose the motions. Two peculiar phenomena, namely, long (perhaps super-persistent) chaotic transients, and fluctuating Lyapunov exponents, are observed; in the two such cases the LLE method fails to work. In addition, the effects of various system parameters, namely, the position of the elastic axis, the frequency ratio, the airfoil/air mass ratio, the viscous damping ratios, and the location of the center of mass, on the response of the aeroelastic system, are investigated.