A numerical method for the variable-order fractional functional differential equations (VO-FFDEs) has been developed. This method is based on approximation with shifted Legendre polynomials. The properties of the latter were stated, first. These properties, together with the shifted Gauss-Legendre nodes were then utilized to reduce the VO-FFDEs into a solution of matrix equation. Sequentially, the error estimation of the proposed method was investigated. The validity and efficiency of our method were examined and verified via numerical examples.