Several comments first:
You should use is, not evalb.

You need more assumptions, like positivity of other constants.

The name gamma stands for Euler's constant which is roughly 0.5772156649.
I would recommend using assuming instead of assume.
This works:

restart;
is((a - c)/(3*b) < gamma*(a - c)/(3*b*gamma - 1)) assuming 2*c < a, 0 < a, a/c < 3*b*gamma,c>0,b>0;
### Answer true

There is no error:
After all int(f(x),x=0..1) is exactly the same as int(f(y),y=0..1).
Thus I see nothing wrong here:

J := Int(r[1]^2*phi[1](r[1]), r[1] = -infinity .. infinity)
     *
     Int(r[2]^2*phi[2](r[2]), r[2] = -infinity .. infinity); 

simplify(J);
##Result: Int(r[1]^2*phi[1](r[1]), r[1] = -infinity .. infinity)*Int(r[1]^2*phi[2](r[1]), r[1] = -infinity .. infinity)
## There is only visual difference between the above result and that from:
map(simplify, J);
## The name of the integration variable in a definite integral is irrelevant.

Here is a much simpler looking example:

restart;
J:=Int(f(x),x=0..1)*Int(g(y),y=0..1);
simplify(J);
simplify(value(J));

Use unevaluation quotes around the whole expression:
 

'P(X <= 5) = P(-5 <= -X) and P(-5 <= -X) = P(E(X) - 5 <= E(X) - X)';

The evaluation to false will come after a % is executed as in
 

'P(X <= 5) = P(-5 <= -X) and P(-5 <= -X) = P(E(X) - 5 <= E(X) - X)'
%;  # false

evalb(2+2=4);
https://mapleprimes.com/questions/237373-How-To-Answer-Your-Own-Question-On-MaplePrimes

This is rather intriguing.
It seems that the difference lies in the different addresses of your tabletop produced array and the one coming out of Maple's DEplot:
 

restart;
ode1:=2*y(t)+t*diff(y(t),t) = t^2-t+1;
p1:=DEtools:-DEplot(ode1,y(t),t=0..3.5,y=0..3):
#####################################
ode2:=3*y(t)+t*diff(y(t),t) = t^2-t+1;
p2:=DEtools:-DEplot(ode2,y(t),t=0..3.5,y=0..3):
#####################################
A0:=Array(1 .. 3,1 .. 2,{(1, 1) = HFloat(undefined), (1, 2) = HFloat(undefined), (2
, 1) = HFloat(undefined), (2, 2) = HFloat(undefined), (3, 1) = HFloat(undefined
), (3, 2) = HFloat(undefined)},datatype = float[8],order = C_order);
#####################################
whattype(A0);
has(A0,HFloat(undefined)) ;
A1:=op([1,2],p1);
whattype(A1);
has(A1,HFloat(undefined)) 
A2:=op([1,2],p2);
whattype(A2);
has(A2,HFloat(undefined))
######################################
H0:=indets(A0,float)[1];
addressof(H0);
H1:=indets(A1,float)[1];
addressof(H1);
H2:=indets(A2,float)[1];
addressof(H2);
addressof(H1)-addressof(H2); # 0
has(A0,H0); #true
has(A1,H1); #true
has(A2,H2); #true
has(A0,H1); #false
has(A2,H1); #true

If you continue in the same session with this:
 

HFloat(1/0.);
addressof(%); 

you get an address different from the previous.

The code:

int( (3*x^2-2*x+1)/sqrt(x^3+x^2-2*x),x);
simplify(%);

The result of the simplification is this:
2*(2*sqrt(-x)*(EllipticE(1/2*sqrt(2*x + 4), 1/3*sqrt(6)) - 1/12*EllipticF(1/2*sqrt(2*x + 4), 1/3*sqrt(6)))*sqrt(-3*x + 3)*sqrt(2)*sqrt(2*x + 4) + x^3 + x^2 - 2*x)/sqrt(x^3 + x^2 - 2*x)

### If you want definite integrals, you do like this:
 

f:=(3*x^2-2*x+1)/sqrt(x^3+x^2-2*x); # For convenience
plot(f,x=-5..5); # Plotting
## Notice that nothing is plotted between 0 and 1.
## That is because f is imaginary in that interval.
int(f,x=0..1); # Exact
evalf(%); #  A numeric approximation
int(f,x=0..1,numeric); # Directly
int(f,x=1..10,numeric);
int(f,x=1..10);
evalf(%);

Since this initial value problem has infinitely many solutions it would cost some more work to write down a formula for these.

The solution we receive first y(x) = 0 is the one that dsolve/numeric would have to come up with.

restart;
ode:=diff(y(x), x) = sqrt(y(x))*sin(x);

dsolve([ode,y(0)=0]);

sol:=dsolve([ode,y(0)=0],'implicit');
eq:=isolate(sol,sqrt(y(x)));
sol2:=op(solve(eq,{y(x)})) assuming y(x)>=0; #sqrt(z) returns the principal branch.
plot(rhs(sol2),x=-2*Pi..2*Pi);
## An example of the infinitely many solutions (not identically zero) to the initial value problem:
sol3:=y(x)=piecewise(x>=0 and x<2*Pi,0,rhs(sol2));
plot(rhs(sol3),x=-2*Pi..4*Pi);
####
res:=dsolve([ode,y(0)=0],numeric);
plots:-odeplot(res,[x,y(x)],-2*Pi..2*Pi,thickness=3);

There is no problem with invlaplace in my Maple 2023.1:
kernelopts(version);
`Maple 2023.1, X86 64 WINDOWS, Jul 07 2023, Build ID 1723669`

You have the code, so what's the problem:
 

restart;

ode:=diff(y(x),x)=1-2*y(x/2)^2;                                                                                                       
ics := y(0) = 2;

y3:= dsolve({ics, ode},numeric, delaymax = 1.7);

YY5:=plots:-odeplot(y3,x=0..10); 

This has been the case forever. Unassigned names are treated as nonzero by many Maple commands.
Take the simple examples:
 

restart;
solve(a*x=0,x); # 0
evalb(a=0); # false
is(a=0); # false
is(a=0) assuming a=0; # true, obviously
1/a;
1/0; # error

For solve there is these days the parametric option:
 

solve(a*x=0,x,parametric);

In your case solve/parametric could be used like this:

A:=Matrix([[1, 0, 0, 0, a__1], [0, 1, 0, 0, b__1], [-216, -18, 0, 0, c__1], [0, 0, 1, 0, d__1], [0, 0, 0, 1, e__1]]);
sys:=LinearAlgebra:-GenerateEquations(A,[x1,x2,x3,x4,x5],<0,0,0,0,0>);
solve(sys,[x1,x2,x3,x4,x5],parametric);

The answer shows correctly that the condition for the zero solution is:
1/216*c__1 + 1/12*b__1 + a__1 <> 0

I would begin by plotting your expression y which is clearly even, but imaginary for 1<x <2.
You don't need any readlib these days. It has been unnecessary for many years.

restart;

y:=sqrt((x^2*(x^2-4))/(x^2-1));
###
plot(y,x=0..10);
eval(y,x=3/2);
evalc(%);
limit(y,x=infinity);
limit(y,x=1,left);
ymax:=maximize(y,x=2..10);
ymin:=minimize(y,x=2..10);
ymin:=minimize(y,x=0..1-1e-10);

The expression  is singular at (x,y)=(0,0):
eval(R0,{x=0,y=0}); ## error
Try another point e.g. (x,y) = (1,0):
 

mtaylor(R0, [x=1, y=0], 5);

I think I managed to read your ode.
I would strongly advocate using odeplot instead of phaseportrait.
You cannot get a direction field (arrows) from odeplot, but your ode is non-autonomous, so a direction field has no meaning anyway.
 

restart;
ode:=diff(y(x),x,x)+8*diff(y(x),x)+9*y(x)+36*y(x)^2=2*sin(x);
## A plot of y(x) versus x:
DEtools:-phaseportrait(ode,y(x),x=-0.2..5,[[y(0)=1,D(y)(0)=0]],numpoints=1000);
# Rewriting ode as a first order system:
ode1:=diff(y(x),x)=v(x);
ode2:=subs(ode1,ode);
## A plot of diff(y(x),x) (i.e. v(x)) versus y(x):
DEtools:-phaseportrait([ode1,ode2],[y(x),v(x)],x=-.2..5,[[y(0)=1,v(0)=0]],numpoints=1000);
#### Much simpler to use dsolve/numeric directly on ode:
res:=dsolve({ode,y(0)=1,D(y)(0)=0},numeric);

## Now plot with odeplot:
plots:-odeplot(res,[y(x),diff(y(x),x)],-.2..5,numpoints=1000,thickness=2);

You may have a look at the help page ?evalf,details.
Also to get inspired you could try showstat(`evalf/constant/Catalan`).
Here is a rather silly example:

restart;
`evalf/constant/K`:=proc()  evalf(sqrt(2)*sin(sqrt(17))) end proc:
evalf(K);
evalf[15](K);
evalf[20](K);
evalf(exp(K)*Pi);
#################
showstat(`evalf/constant/Catalan`);

Since I believe that Maple's solve simply uses a formula for solving polynomials of degree less than 5 the answer that comes out from solve will just be the result of inputting the coefficients of the polynomial into the formula.
Using the very same polynomial as mmcdara, here is a graphical illustration, where the returned sequence from solve is kept in a list, not a set. Thereby we are not relying on any special knowledge of the ordering of sets in Maple.
In the grapical illustration you see that none of the 3 solutions remain real: The first in the list is red, the second blue, and the third is green.

restart;
pol:=x^3+a*x+1;
sol:=[solve](pol,x); # result: a list
with(plots):
complexplot(sol,a=-20..20,color=[red,blue,"DarkGreen"],thickness=3,style=point);
animate(complexplot,[sol,a=-20..A,color=[red,blue,"DarkGreen"],thickness=3,style=point], A=-20..20);

The animation follows   :