The Blume-Emery-Griffiths model is simulated using the cooling algorithm which is improved from the Creutz cellular automaton (CCA) under periodic boundary conditions. The simulations are carried out on a simple cubic lattice at K/J = -1.5 in the range of -3.5 < D/J < 0.5, with J and K representing the nearestneighbour bilinear and biquadratic interactions, D being the single-ion anisotropy parameter. The phase diagram characterizing phase transition of the model is obtained. We found different kinds of phase transitions between the ferromagnetic, quadrupolar, staggered quadrupolar and ferrimagnetic phases for K/J = -1.5. In particular, the region of the phase diagram containing a ferrimagnetic phase is explored and compared to those obtained by other methods. The simulations confirm that the ferrimagnetic phase occurs in the narrow interval -3.006 <= D/J < -3. This result is in a good agreement with Monte Carlo renormalization group and closer to the cluster variation method result than the mean field approximation result.