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).

OK, I have installed the update and everything seems to be OK, no problems detected (Win7, 64 bit).
The mentioned cmaple help problem is still present but it's not a very serious bug.

About your 2D example in the help page issue: the argument of value should be a label; you have probably unchecked "Show equation labels" In Tools > Options > Display. But even so, value(...) executes correctly after pressing Enter.
Switching 1D <--> 2D also works.

It is true that 1D <--> 2D  does nothing in the "How do I" page, but note that this is not actually a help page, but an embedded document. It seems that for all real help pages it works.

Maple has very good implementations for both, see  ?SmithForm,  ?HermiteForm

They even work in general Euclidean domains, see  ?LinearAlgebra[Generic]:-SmithForm

You just need:

Eigenvectors(K, evalf(M));

(omit evalf for a symbolic solution).