You could look at this:

restart;
showstat(`convert/rational`);
interface(verboseproc=2);
eval(`convert/rational`);

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`);