The temperature dependence of the magnetization and quadrupole order parameters of the Blume-Emery-Griffiths (BEG) model Hamiltonian with the nearest-neighbor ferromagnetic exchange interactions [both bilinear (J) and biquadratic (K)] and crystal field interaction (D) is studied using the lowest approximation of the cluster variation method. Besides the stable solutions, metastable and unstable solutions of the order parameters are found for various values of the two different coupling parameters, alpha = J/K and gamma = D/K. These solutions are classified using the free energy surfaces in the form of a contour map. The phase transitions of the stable, metastable and unstable branches of the order parameters are investigated extensively. The critical temperatures in the case of a second-order phase transition are obtained for different values of alpha and gamma calculated by the Hessian determinant. The first-order phase transition temperatures are found using the free energy values while increasing and decreasing the temperature. The temperature when both the free energies equal each other is the first-order phase transition temperature. Finally, the results are also discussed for the Blume-Capel model which is the special case of the BEG model. (C) 1999 Elsevier Science B.V. All rights reserved.