Wrinkles in layered neo-Hookean structures were recently formulated as a Hamiltonian system by taking the thickness direction as a pseudo-time variable. This enabled an efficient and accurate numerical method to solve the eigenvalue problem for onset wrinkles. Here, we show that wrinkles in graded elastic layers can also be described as a time-varying Hamiltonian system. The connection between wrinkles and the Hamiltonian system is established through an energy method. Within the Hamiltonian framework, the eigenvalue problem of predicting wrinkles is defined by a series of ordinary differential equations with varying coefficients. By modifying the boundary conditions at the top surface, the eigenvalue problem can be efficiently and accurately solved with numerical solvers of boundary value problems. We demonstrated the accuracy of the symplectic analysis by comparing the theoretically predicted displacement eigenfunctions, critical strains, and wavelengths of wrinkles in two typical graded structures with finite element simulations.