Summary: This work proposes a mathematical model and its parallel implementation on two different computer architectures to simulate three-dimensional bursting phenomena. The mathematical model consists of four, nonlinearly coupled partial differential equations and includes fast and slow subsystems. The differential equations have been discretized by means of a linearly-implicit finite difference method in equally-spaced grids. The resulting system of linear algebraic equations at each time level has been solved by means of the Generalized Minimal Residual (GMRES) method with Jacobi preconditioning. The parallel implementation has been formulated using the message passing paradigm, where the linear system of equations has been solved by means of the parallel sparse solver GMRES included in the portable extensible toolkit for scientific computation (PETSc). The parallel simulation has been evaluated on (1) A cluster of biprocessors Xeon (TM) and (2) A SGI Altix 3700 Bx2. Similar performance has been obtained on both platforms, with better scalability on the DSM multiprocessors.