This study presents a shifted Bernstein polynomial-based method for numerically solving the variable fractional order control equation governing a viscoelastic bar. Initially, employing a variable order fractional constitutive relation alongside the equation of motion, the control equation for the viscoelastic bar is derived. Shifted Bernstein polynomials serve as basis functions for approximating the bar’s displacement function, and the variable fractional derivative operator matrix is developed. Subsequently, the displacement control equation of the viscoelastic bar is transformed into the form of a matrix product. Substituting differential operators into the control equations, the control equations are discretized into algebraic equations by the method of matching points, which in turn allows the numerical solution of the displacement of the variable fractional viscoelastic bar control equation to be solved directly in the time domain. In addition, a convergence analysis is performed. Finally, algorithm precision and efficacy are confirmed via computation.