So we talked about this on the list a while ago (see http://groups.google.com/group/sympy/browse_thread/thread/40cafb75db3cb522/eff6ff30dcd450f9? lnk=gst&q=RootOf). The idea is to have solve() return some kind of RootOf class whenever it cannot
solve an equation instead or raising NotImplemented or returning [] as it does now. Both Maple and
Mathematica do this. Maple uses RootOf and Mathematica uses Solve. I think we should use RootOf,
because we already have that class for Polys.

The RootOf class should probably be pulled out of the Polys module, as we will want to use it for all
equations, not just polynomial ones. Here is the link to the Maple docs on how they do it: http://www.maplesoft.com/support/help/view.aspx?path=RootOf/indexed and http://www.maplesoft.com/support/help/view.aspx?sid=19556 . I don't own Mathematica and can't
find the link to their docs on Solve(). Does anyone know it?

In Maple anyway, there are two kinds of RootOf's, indexed and unindexed. Indexed RootOf's are
solutions to polynomials, where you know how many roots there are. These are what the current
RootOf's in SymPy are, but presently, they are not well defined. The indices should be defined in a
manner that you can will always get the same root for the same index with evalf. Here is from the
Maple docs (see link above) on how they order the roots:

1. The equation expr is a polynomial in x whose coefficients can be numerically evaluated (see
   complexcons). The RootOf represents the ith complex root (multiplicities taken into account) of the
   polynomial. The roots are ordered counter-clockwise, with ties broken by increasing modulus. The first
   root is chosen as follows:

• roots with the smallest argument in absolute value are selected,
• roots with positive argument, if any, are then chosen,
• the root with the smallest modulus amongst the remaining roots is the root number 1.
  If the RootOf represents an nth root, that is, if expr is a binomial of degree n, then the index i
  corresponds to exp((i-1)2I*Pi/n)*expr^(1/n) where expr^(1/n) is the principal branch of the nth root
  function.
  The default ordering may be changed by over-writing the function called RootOf/sort. The argument
  of this function is a list of complexcons and the result must be the list of sorted elements. If the default
  ordering is modified, some functions may not work as expected (convert,radical for instance). The evala
  package will still behave as expected.

I think we should use whatever order is most efficient for evalf. The ability to evalf() RootOf's is also
important. They would be quite useless without this.

There should also be a way to convert a RootOf into an equality. So if you have RootOf(cos(x) - x, x),
you could convert it into Eq(cos(x) - x, 0) or Eq(cos(x), x). There should be a way to do this with nested
RootOfs, at least for polynomial RootOfs (Maple can do it anyway).

For non-indexed RootOf, I think we should use RootsOf as Ondrej suggested on the list. I don't think
to use the same ordering as indexed roots (consider it the expression had in infinite number of roots
that were all over the complex plane). Maple lets you use a bounding box from two complex number,
or you can label RootOfs.

We should make re() and im() work, so that im(RootOf) returns 0 if we know that root is strictly real.

Also, we will need to fix issue 4656 .

Also, I wonder how this would word for systems of equations. I would like for it to work for that too.
For example, if you try to solve for the arbitrary constants in an ODE with initial conditions and solve()
can't do it, then you could return RootOf's and the user could then evalf if they wanted to.