Typically galaxies contain 1011 stars and globular clusters contain 106 stars. In order to investigate the dynamical evolution of these systems one may perform numerical simulations (N-body simulations); which in turn are restricted to 103 to 107 particles depending on the required accuracy and adopted evolutionary model. For the N-body direct method the computational effort scales at least as O(N3), and thus the need for efficient and accurate evolutionary models and advances in hardware in apparent. A review of the N-body evolutionary models is presented and, with the advent of parallel super computers, the parallelisation of these evolutionary models is investigated. In particular, we discuss the parallelisation of a direct method N-body code; focusing upon a new parallelisation strategy, implementing the evolutionary model with a portable parallel language (High Performance Fortran), and reducing the communication costs between processors by a suitable implementation of the Hypersystolic algorithm. A dynamical system, such as a galaxy, may oscillate about a stable equilibrium if it is excited above that state. These modes of oscillation may persist long enough to have observable consequences, despite being weakly damped. The moods of oscillation are investigated with a range of evolutionary models. The simplest mode to detect and simulate is the "fundamental" mode or "l = 0" mode, which manifests itself as a radical oscillation of the entire system. To investigate this mode a King W0 = 1 theoretical model is adopted. The effect of adopting a softened gravitational potential to generate the forces between stars on the fundamental mode is investigated with a direct method N-body evolutionary model, and these results are compared to those which use an SCF (self-consistent field) type code. The presence of small and amplitude fundamental mode oscillations is detectable when the perturbation particle method is used, which would otherwise be undetectable for another N-body evolutionary model, due to particle noise. Furthermore, the source of heavy damping in direct method N-body simulations is found to be phase mixing. The experimental fundamental mode oscillations are found to match well with the theoretically predicted frequencies. A more complicated mode, the "sloshing" mode or "l = 1" is also studied. This mode manifests itself in the density centre shifting or sloshing about. To investigate this mode a King W0 = 5 model is adopted, and the results compared to an analytical predicted frequency.