It is well known that an Earth-like planet has four stationary solutions, two stable and the other two unstable in the linear sense, which can be checked by computing the eigenvalues of the linearized system of equations of motion. Determining the orbital stability (or Liapunov stability) is rather more involved, and it requires firstly, the computing of the normal form around the equilibrium; secondly the expression of this normal form in action-and-angle variables in order to apply the so called Arnold's theorem of stability. For some specific values, higher order normal forms are required. However, Arnold's theorem is useless in the presence of resonances, and a new technique (related again with normal forms) must be applied in order to determine orbital stability. In this work, we proceed symbolically, taking the harmonic coefficients as parameters. We find the Liapunov stability diagram in the parametric plane for the stationary points. Resonances 2:1 and 3:1 are also studied, as well as the case at which a higher order normalization is required. The advantage of the symbolic analysis is that it is enough to replace the actual values of a planet or celestial body to have the orbital stability of the stationary points.