@nm Mathematica is not the supreme authority.

P.S. With your statement,  (x + y(x))^7 * (1+diff(y(x),x)) = 0  has 8 solutions.

@nm It's not in the standard form (or reducible to)  d/dx y(x) = F(x, y(x)).

The reason is obvious. Anyway, I see ODESteps as designed for standard simple ODEs. 

@nm No bug without Physics updates.

@minhthien2016 Insert f:=simplify(f) assuming real; before calling extrema. 

@Carl Love Very nice, vote up! Unfortunately the rose is not quite 3d: rotating the plot, several gaps are visible and the stem is just "separately added".

@Hullzie16  Adding the option numeric, dsolve produces a nrmerical solution,
Your system has a lot of floats, so, approximations. In such situations, an exact (analytic) solution is usually useless (even when it exists). 

@Ali Hassani 

seriestorec(s, A(n)) 
   ==> for a given series s, finds a "suspected" recurrence relation of its coefficients A(n)
rectodiffeq(rec, A(n), f(z))
   ==> for a recurrence rec, finds the ode satisfied by the "generating function"
           f(z) = Sum(A(n)*z^n, n = 0 .. infinity).
 

@dharr  When checking the minimal polynomial we must use a much higher precision.

Digits:=200:
evalf(expr) - fsolve(poly)[2]: evalf[15](%);
#                        1.37954446875327 e-38

(it should be near 1e-200).
So, the first poly is incorrect.

The second one (of degree 20) is ok.

@sursumCorda  Why not use e.g.

Matrix( [ [x || (1 .. 2)] , [y || (1 .. 2)] , [1 $ 2] ] );

Download W.mw

@nm Maple decided to use the -1 branch because it knows the LambertW function. For example,  the restrictions  f₁:(-∞,-1]→[-1/e,0)  and  f₂: [-1,∞)→[-1/e,∞)  of  f(x)=x*exp(x)  are inverted as LambertW(-1,y) and LambertW(0,y). Given IC, Maple knows which one to use.

The same thing happens in other situations, e.g to invert  cos : [4*Pi, 5*Pi] -> [-1,1]. So it's nonsense to look for an algorithm, because it would imply to rewrite Maple.

If you do not want IC, just use the option implicit as I said; most cases will be solved this way.

@C_R Actually, in both situations, solve only finds a few surfaces contained in the interior of the torus (the main components of the solution are missing).

@Preben Alsholm The result of

sol2:=SolveTools:-SemiAlgebraic({torus1<0}, [x,y,z]);

seems to be correct but more complicated than necessary.

After completing my last post I was said that there are no tags; but they were introduced! I accepted a proposed tag (a stupid one) and the post was accepted.