E.g. (by hand):

It seems that you are using plotsetup(char);

To restore the defaults, use

plotsetup(default);

You must be aware that a series structure is obtained only using the series command, e.g.

series(exp(z),z,3);
    1+z+(1/2)*z^2+O(z^3);

But if your input is

1+z+(1/2)*z^2+O(z^3);

this is no longer a series; it's a sum of 4 terms and O is now simply a symbol.

Edit. Note also that the command  MultiSeries:-multiseries  uses another structure named SERIES; this is more general because it can use other "scales".

fsolve finds multiple roots only for polynomials.
In your case you may use solve (because the equation can be solved exactly i.e. symbolically) or RootFinding
(there are other possibilities).

Download rootfinding.mw

Try to include the path in the filename.

P:=sort(E,output=permutation);
E__sorted:=E[P];
Q__sorted:=(Q^+)[P]^+;

   Use   dsolve.

You must first express it in terms of y1 and y1'.

y1tt := rhs(isolate(eq1, diff(y[1], t, t))):
plots:-odeplot(sol,[t,y1tt],t=0..10);

convert inserts \n  and it is correct from Maple's point of view because white spaces are ignored by the parser.

You can remove them saving e.g. the string
StringTools:-DeleteSpace(result_as_string):

1) procedures use last name evaluation, so, eval must be used.
2) nonIdMaps is [f,f]  and now  f is the last index of the loop, i.e.  x -> 3*x
 

nonIdMaps := []:
for f in [x -> x,x -> 2*x,x -> 3*x] do
   eval(f),f(y);
   if f(y) <> y then nonIdMaps := [nonIdMaps[],eval(f)] end if
end do;
nonIdMaps,map(f -> f(y),nonIdMaps);

works as expected.

Your t should be x (probably, otherwise x is a parameter).
Do you have an initial condition? It seems that the ODE has not real solutions (y should be >0).

I very seldom use Physics, but it seems that a prefix acts only if the rest of the name is  n  or  __n,
n containing only decimal digits.

Obviously y > 0.

a) For x>0 the inequality reduces to some

x^2 < f(y).

b) For x<0 the signum of the denominator   x+0.015*y^(1.2)  must be considered too.

It should be easy to start from here.

I think that RealDomain has a very limited applicability. I prefer not to use it.
It redefines several standard functions. For example,

exp(I*t);
       undefined

Why do you want to use it with Statistics?
 

[1]. Your problem does not seem to be a "pure" mathematical one; after all you have "floats".
So, a global minimization/maximization package (such as DirectSearch) should be enough.

[2]. Suppose now that you still want to check the concavity.
I would check that the Hessian is <= 0 for the variables in a "dense" subset of the domain (this should be enough for a plausible conclusion). It is also possible to take a consistent set of randomly generated points.
Unfortunately at least the first function (sol1) is not concave:

H1:=eval(Hess_TP1, {E=1, T=0.2, W=2, p=21} ):
IsDefinite(H1, 'query' = 'negative_semidefinite');
        false

You could use these suggestions for other parameters if you want, but I'd recommend [1] (anyway the cancavity is not necessary in many cases).