Many physical models require the simulation of a large number ($N$) of particles interacting through pair-wise inverse square law forces. $N$-body simulations are employed in fluid-dynamics, biochemistry, astrophysics, electrodynamics and molecular dynamics. The computational problem is intrinsically hard and these simulations are time-intensive. Existing algorithms exploit either the spatial proximity of particles or the temporal proximity of states. In this thesis, we formally combine the two approaches and present an algorithm with sequential time complexity $O(N^{4/3})$ to integrate $N$ uniformly distributed particles in 3D over one crossing time against the $O(N^{8/3})$ complexity of the direct method. Under reasonable assumptions, our algorithm is optimal. The core of the algorithm is the temporal multipole expansion of the field in terms of the space-time coordinates of the field-point. We also present efficient parallel algorithms on the 2D and 3D Mesh and Hypercube which amortize communication costs through temporal multipole expansions. The parallel algorithms offer an order of magnitude improvement over existing algorithms for even $10^{4}$ particles. A sequential implementation of the algorithm for two-dimensional $N$-body systems shows the predicted asymptotic scaling. A parallel version on a 16-processor Intel iPSC/860 machine is also in conformance with theoretical expectations.