We consider fast computation methods for simulation of dynamics of a collection of particles dispersed in an unbounded Stokesian suspension. Stokesian suspensions are of great practical interest in the manufacturing and processing of various commercial products. The most popular dynamic simulation method for these kind of suspensions was developed by Brady and Bossis (Brady and Bossis [1988]). This method uses a truncated multipole expansion to represent the fluid traction on particle surfaces. The hydrodynamic interactions in Stoke-sian suspension are long ranged in nature, resulting in strong coupled motion of all particles. For an N particle system, this method imposes an O(N3) computational cost, thus posing limitations to the number of particles that may be simulated. More recent methods (Sierou and Brady [2001], Scintilla, Darve and Shaqfeh [2005]) have attempted to solve this problem using Particle Mesh Ewald summation techniques by distributing the moments on a grid and using Fast Fourier Transform algorithms, resulting in an O(N log N) computational cost. We review these methods and propose a version that we believe is some-what superior. In the course of this study, we have identified and corrected errors in previous studies that maybe of some importance in determining the bulk properties of suspensions. Finally, we show the utility of our method in determining certain properties of suspensions and compare them to existing analytical results for the same.