Not tremendously difficult, but not nearly as simple as assuming circular orbits. The mechanics book by Kleppner and Kolenkow has a good derivation of this.
They can also be parabolic or hyperbolic, if the planet has a sufficient escape velocity. Finding the allowed orbits from Lagrange's equations of motion may require some work, unless you already know the solution (in which case you can test it by substitution).
I imagined the following approach. The equation of an ellipse is ##x^2/a^2 + y^2/b^2 = 1##. If we define ##\bar{x} = x / a, \bar{y} = y / b## we have ##\bar{x}^2+ \bar{y}^2 = 1## which is the equation of a unit sphere.
Of course that satisfy Newton's equation of gravity and so do ##x^2 / a^2## and ##y^2 / b^2##.