The Legendre polynomial series approach (LPSA) has been widely used to solve guided wave propagation in various structures since 1999. The LPSA directly introduces the boundary conditions into the control equations through the rectangular window function. However, in the solving process, the Legendre polynomial series, the rectangular window function, and their derivatives are introduced into the integral kernel functions, which results in a lot of CPU time on the abundant numerical integration calculations. To overcome this defect, an analytical integration Legendre polynomial series approach (AILPSA) is presented and is used to solve the Lamb waves in fractional order thermoelastic multilayered plates. Coupled wave equations and heat conduction equation are solved by the AILPSA and the LPSA, respectively. Comparison between two approaches indicates the computational efficiency of the AILPSA is improved by more than 90%. In addition, a backward error estimation is given in the convergency analysis to make up for the deficiency of the numerical experiments. Finally, the influence of fractional order on the thermoelastic Lamb wave is discussed.