In this paper, an isothermal compositional grading process is simulated based on convex splitting methods with the Peng-Robinson equation of state. We first present a new form of gravity/chemical equilibrium condition by minimizing the total energy which consists of Helmholtz free energy and gravitational potential energy, and incorporating Lagrange multipliers for mass conservation. The time-independent equilibrium equations are transformed into a system of transient equations as our solution strategy. It is proved our time-marching scheme is unconditionally energy stable by the semi-implicit convex splitting method in which the convex part of Helmholtz free energy and its derivative are treated implicitly and the concave parts are treated explicitly. With relaxation factor controlling Newton iteration, our method is able to converge to a solution with satisfactory accuracy if a good initial estimate of mole compositions is provided. More importantly, it helps us automatically split the unstable single phase into two phases, determine the existence of gas-oil contact (GOC) and locate its position if GOC does exist. A number of numerical examples are presented to show the performance of our method.