Question: Maximize expression with RootOf and parameters

I have equations of the form x=max( f(x,u,E) ), where x is the unknown "variable" to solve for, where u is the "control" variable that can be selected to maximize f, and where E is a parameter. These are all real, positive scalars on [0,1].

One of my objectives is to plot x against E, for those values of u that solve the maximization. In my problem, the optimal values of u are expected to be unique for a given value of E.

My approach was to define a grid of values for the parameter E and to maximize f(x,u) subject to the constraint that x=f(x,u), for each value of E in the grid (here I dropped E as an argument of f to indicate that E has been assigned a specific value). 

For that purpose, I used Maximize from the Optimization package. I could make it to work on simple cases, but ran into problems for less simple cases. 

It will almost certainly be clearer if I give examples.

Example 1: Simple case: works
restart;
eq := x = (.6*E^4*u^4+.4*E^4*u^5-8000.*x^5)/(E^4*u^4);
xsol := solve(eq,x);
seq(Optimization:-Maximize(xsol,u=0..1,'variables'=[u],assume=nonnegative),E=0..1,0.1);

Example 2: Less simple: fails
restart;
eq := x = 11/500+6*u*(1/25)-u^2+(-174200000000*x^5+13045131*E^4*u^4)^(1/5)/(65225655*E^4*u^4);
xsol := solve(eq,x) assuming x>0 and x<1 and u<1 and u>0 and x::positive:
seq(Optimization:-Maximize(xsol,u=0..1,'variables'=[u],assume=nonnegative),E=0..1,0.1);
Error, (in content/polynom) general case of floats not handled

I'm not sure what this message means. I wrapped the equation with evalf() in case the problem was related to the fact that it's specified in terms of rationals rather than floats, but that didn't help.

If I fix the parameter E to a single value instead of a list of values, I get another error message:
Optimization:-Maximize(eval(xsol,E=1),u=0..1,'variables'=[u],assume=nonnegative);
Error, (in Optimization:-NLPSolve) complex value encountered

However, a quick plot shows that there is a maximum (obviously) and that the curve is not flat but quite strongly curved, suggesting to me that it should be easier for Maximize to find the max.

Once I know where the maximum is located, I can change the range of values to search over and finally find the maximum: 

Optimization:-Maximize(eval(xsol,E=1),u=0.1..0.3,'variables'=[u],assume=nonnegative);
[0.0360041440617226, [u = 0.126172555344267]]

But I'd like a systematic approach that does not rely on tweaking, because I would like to run the simulation multiple times on different kinds of functions. 

One idea I had was to use try/catch to grope for a range for E, starting with 0..1, then upping it to 0.1..1 in case an error was "caught", etc. Before I get started on this "project" I'd like to hear opinions! There may be a more straightforward approach.

Any suggestions will be much appreciated.


plot( eval(xsol,E=1),u=0..1);