A novel recursive aggregation algorithm for numerical simulations of dynamic systems is proposed and analyzed. The algorithm exploits a special structure of the linear equation problem resulting from the discretization of the dynamic system and an aggregation/disaggregation procedure. The algorithm has a time complexity of (I(q) + 2M(q) + 3) log N for solving linear systems with q states for N discrete time instants, using O(q3 N) processors, where I(q) is the parallel time complexity for inverting a q x q matrix, M(q) is the parallel time complexity for matrix multiplication of two q x q matrices. The competing parallel cyclic reduction method for the same problem has a time complexity of (I(q) + 3M(q) + 4) log N. Thus, the proposed algorithm has a definite speed advantage over the cyclic reduction method. An approximation technique for the unknown boundary conditions in boundary value problems is also proposed. The algorithm was implemented to simulate some dynamical (stable and unstable) systems, and the numerical results show that the accumulation of roundoff errors is insignificant as compared to the discretization errors.