Equilibrium properties of BlumeâEmeryâGriffiths (BEG) model with the arbitrary bilinear (J), biquadratic (K) and crystal field interaction (D) are studied on a body centered cubic lattice by using the pair approximation of the cluster variation method, which is identical to the Bethe approximation. The thermal variation of the stable, metastable and unstable dipole and quadrupolar moment order parameters are found for various values of the two different coupling parameters J/K and D/K. 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 J/K and D/K by the Hessian determinant. The first-order phase transition temperature is determined by matching the values of the two branches of the free energy values. On the other hand, nonequilibrium behavior of the model is studied by using the path probability method with the pair distribution and the set of nonlinear differential equations, which is also called the dynamic or rate equations, is obtained. The solution of the dynamic equations is expressed by means of flow diagrams. The stable, metastable and unstable solutions are shown and the "overshooting" phenomenon is seen in the flow diagrams, explicitly. The role of the unstable points, as separators between the stable and the metastable points, is described and how a system freezes in a metastable state is also investigated, extensively. Moreover, the dynamic equations are also solved by using the RungeâKutta method in order to study the relaxation of order parameters. In this investigation, the "flatness" property of metastable states and the "overshooting" phenomenon are seen explicitly.