Aunque el cálculo fraccionario ha evolucionado significativamente desde su origen en la correspondencia de 1695 entre Leibniz y L"Hôpital, el tratamiento numérico de ecuaciones diferenciales fraccionarias de orden múltiple sigue siendo un desafío. Los métodos existentes son a menudo costosos computacionalmente o dependen de marcos operativos complejos que obstaculizan su aplicabilidad más amplia. En el presente estudio, se propone un nuevo algoritmo numérico basado en funciones híbridas ortogonales (HFs), que se construyeron como combinaciones lineales de funciones de muestreo y retención constantes por tramos y funciones triangulares lineales por tramos. Estas funciones, que pertenecen a la clase de polinomios ortogonales de grado 1, se emplearon para obtener la solución numérica de ecuaciones diferenciales fraccionarias de orden múltiple definidas en el sentido de Caputo. Se derivó una matriz operativa generalizada de un solo paso para expresar explícitamente la integral fraccionaria de Riemann-Liouville de las HFs en términos de las propias HFs. Esto permitió transformar la ecuación diferencial fraccionaria de orden múltiple original directamente en un sistema de ecuaciones algebraicas, simplificando así el proceso de solución. El algoritmo desarrollado se aplicó a una serie de problemas de referencia, incluyendo tanto FDEs lineales como no lineales de orden múltiple con coeficientes constantes y variables. Las comparaciones numéricas con métodos bien establecidos en la literatura revelaron que el enfoque propuesto no solo logró una mayor precisión, sino que también redujo significativamente el esfuerzo computacional, demostrando su potencial como una herramienta numérica confiable y eficiente para el modelado de órdenes fraccionarios.