Increasing attention has been given to the application of variable-order fractional differential equations (VO-FDEs) in modeling complex systems exhibiting time-dependent memory and hereditary effects. However, due to the fractional operators being nonlocal and the differentiation order being variable in nature, it remains mathematically difficult to numerically treat boundary value problems associated with VO-FDEs. This paper provides a numerical framework based on symmetry-driven classification for solving a certain class of boundary value problems related to VO-FDEs using a completely symmetric classification procedure via differential characteristic sequence analysis. By utilizing extended symmetries, the original boundary value problem is systematically reduced to an equivalent initial value problem involving ordinary differential equations. This initial value problem can then be solved using a Legendre polynomial operator matrix approach that converts the fractional differential equation into a structured system of algebraic equations via discretization. The results obtained from numerical experiments show that our new approach produces more accurate and faster calculations than the previous methods of solving dynamic systems with fractional derivatives. The use of both symmetry analysis and orthogonal polynomial techniques together to solve these models is confirmed to be a good way to create variable order fractional derivative mathematical models.