In this work we address the dynamic simulation and optimization of chemical processing systems modeled in terms of fractional order differential equations. While fractional derivatives were first proposed by Liouville in 1832 [Samko et al. Fractional Integrals and Derivatives Theory and Applications; Gordon and Breach: New York, 1993; Oldham and Spanier. The fractional Calculus; Academic Press: New York, 1974], fractional differential equation (FDE) models have been only recently been explored. These have been proposed for a wide range of applications that include systems with nonlocal diffusion phenomena and geometries with fractal dimensions. FDE models have been shown to have advantages over traditional integer order models, as they often avoid scale dependent model parameters. For medium or large scale applications of FDEs normally no analytical solutions are available, and therefore, approximated numerical solutions ought to be sought. Moreover, little work has been done to solve fractional order differential equations numerically; most of the existing numerical methods are intended for small scale systems. In this work, we propose a new numerical method, based on Gaussian quadrature on finite elements, to address larger-scale FDEs and extend them to dynamic optimization. To gain a better appreciation about the performance of the new algorithm, we compare its response to recently proposed predictor-corrector methods. We also develop a proposed method to deal with dynamic optimization fractional order systems, an open research problem that has not received wide consideration. We test the performance of the algorithms by deploying three systems embedded with different fractional derivative behavior, from a simple linear dynamic system to a dynamic multiple steady-states bioreactor with various levels of imperfect mixing. The results indicate better numerical properties of the fractional order Gauss quadrature algorithm both for the dynamic simulation and optimal control issues.