In classical study on viscoelastic fluids, researchers commonly ignored the effects of nonlinear convection and dealt only with the cases when the governing equations are linear. This paper presents a research on unsteady boundary layer natural convection heat transfer of Maxwell viscoelastic fluid over a vertical plate. The fractional boundary layer governing equations are firstly formulated and derived. From such derivation, the model constitutes nonlinear coupled equations with mixed time-space derivatives in the convection terms, which are solved by a newly developed finite difference method combined with an L1-algorithm. Results show that both the velocity and the temperature boundary layer manifest a short-term memory and basic relaxation time characteristics. For example, for specified velocity and temperature fractional derivative alpha and beta, the velocity profile has a maximum, which increases and goes from a far-field region to the vertical plate as a decreasing function of a, implying a loss in velocity boundary layer thickness. But the values of such maxima behave conversely for beta. More effects of the fractional order parameters alpha and beta on velocity and temperature fields are shown graphically and analyzed in detail. (C) 2016 Elsevier Ltd. All rights reserved.