In this paper, a kinetic equation of Euler–Bernoulli beam is established with variable order fractional viscoelastic model. An effective numerical algorithm is proposed. This method uses a combination of shifted Bernstein polynomial and Legendre polynomial to approximate the numerical solution. The effectiveness of the algorithm is tested and verified by mathematical examples. The dynamic behavior of viscoelastic beams made of two materials under various loading conditions is studied. [ABSTRACT FROM AUTHOR]