A state space method is proposed for analyzing surface instability of elastic layers with elastic properties varying in the thickness direction. By assuming linear elasticity with nonlinear kinematics, the governing equations for the incremental stress field from a fundamental state are derived for arbitrarily graded elastic layers subject to plane-strain compression, which lead to an eigenvalue problem. By discretizing the elastic properties into piecewise constant functions with homogeneous sublayers, a state space method is developed to solve the eigenvalue problem and predict the critical condition for onset of surface instability. Results are presented for homogeneous layers, bilayers, and continuously graded elastic layers. The state space solutions for elastic bilayers are in close agreement with the analytical solution for thin film wrinkling within the limit of linear elasticity. Numerical solutions for continuously graded elastic layers are compared to finite element results in a previous study (Lee et al., 2008, J. Mech. Phys. Solids, 56, pp. 858-868). As a semi-analytical approach, the state space method is computationally efficient for graded elastic layers, especially for laminated multilayers.