In this paper, we establish a connection between the recently developed data-driven time-frequency analysis \cite{HS11,HS13-1} and the classical second order differential equations. The main idea of the data-driven time-frequency analysis is to decompose a multiscale signal into a sparsest collection of Intrinsic Mode Functions (IMFs) over the largest possible dictionary via nonlinear optimization. These IMFs are of the form $a(t) \cos(θ(t))$ where the amplitude $a(t)$ is positive and slowly varying. The non-decreasing phase function $θ(t)$ is determined by the data and in general depends on the signal in a nonlinear fashion. One of the main results of this paper is that we show that each IMF can be associated with a solution of a second order ordinary differential equation of the form $\ddot{x}+p(x,t)\dot{x}+q(x,t)=0$. Further, we propose a localized variational formulation for this problem and develop an effective $l^1$-based optimization method to recover $p(x,t)$ and $q(x,t)$ by looking for a sparse representation of $p$ and $q$ in terms of the polynomial basis. Depending on the form of nonlinearity in $p(x,t)$ and $q(x,t)$, we can define the degree of nonlinearity for the associated IMF. %and the corresponding coefficients for the associated highest order nonlinear terms. This generalizes a concept recently introduced by Prof. N. E. Huang et al. \cite{Huang11}. Numerical examples will be provided to illustrate the robustness and stability of the proposed method for data with or without noise. This manuscript should be considered as a proof of concept.