This produces the same as your code:
 

restart;
X := Matrix(
   [[1, -X3_3/2 - 1/2, 0, -X2_3], [-X3_3/2 - 1/2, -2*X3_4 - 1, X2_3, 0], [0, X2_3, X3_3, X3_4], [-X2_3, 0, X3_4, 1]]
            );
####
vars := [X3_4, X3_3, X2_3];
w := A^3 - A;

rts:=solve(w,A);
Pols := [(-A^2 + 1)/(3*A^2 - 1), (-A^2 - 1)/(3*A^2 - 1), A*(3*A^2 - 1)*1/(3*A^2 - 1)];
####
vals:=map2(eval,Pols,A=~{rts});
ecs:=[seq](vars=~(vals[i]),i=1..3);
Nelem := [seq(k, k = 1 .. numelems(vals))];
####
Xs:=[seq(eval(X, ecs[i]),i=Nelem)];
####
with(LinearAlgebra):

Eigs := Eigenvalues~(Xs);

BoolEigs:=map((y->is(Im(y) = 0 and 0 <= Re(evalf(y))))~ ,Eigs);

Maple can do it for you. It uses a finite difference method.
 

restart;
eq1 := (diff(f(x), x, x, x))*(a*beta*f(x)^2-1)+(diff(f(x), x))^2-2*a*beta*f(x)*(diff(f(x), x))*(diff(f(x), x, x))+(diff(f(x), x))*(M+k[1])-(diff(f(x), x, x))*f(x)-(alpha*theta(x)+delta*phi(x))/rho = 0;

eq2 := -(diff(theta(x), x, x))*K[SB]*(Df-(Rd+k[hnf]/k[bf])/Pr)+N[t]*K[SB]*(diff(theta(x), x))^2-N[b]*(diff(theta(x), x))*(diff(phi(x), x))-(diff(f(x), x))*(diff(theta(x), x))-lambda*theta(x)-mu*Ec*(M*(diff(f(x), x))^2+(diff(f(x), x, x))^2) = 0;

eq3 := diff(phi(x), x, x)+Le*Sr*(diff(theta(x), x, x))+Le*f(x)*(diff(phi(x), x)) = 0;

ics := f(0) = 0, (D(f))(0) = 0, theta(0) = 1, phi(0) = 1;

bcs := (D(f))(100) = 0, theta(100) = 0, phi(100) = 0;


Parameters1 := rho = 2063.905, k[hnf] = .29942, k[bf] = .2520, mu = .38694, a = .1, beta = 5, k[1] = 2.0, M = 10, alpha = 20, delta = 20, K[SB] = .5, Df = 3, Pr = 1.2, Rd = 5, N[t] = 1.2, N[b] = 1.0, lambda = 1.5, Ec = 5, Le = .1, Sr = .1;

 

sys:={eq1,eq2,eq3,ics,bcs};
res:=dsolve(eval(sys,{Parameters1}),numeric);
with(plots):
odeplot(res,[x,f(x)]);
scene:=[  [x,f(x)],[x,theta(x)], [x,phi(x)]  ];
display(Array([seq(odeplot(res,scene[i],thickness=2,color=blue),i=1..3)]));
scenediff:=evalindets(scene,function,fu->diff(fu,x));
display(Array([seq(odeplot(res,scenediff[i],thickness=2,color=red),i=1..3)]));

You may want a better look at the first part of the interval x = 0..100.
Just add the interval as done here for x = 0..20:
 

display(Array([seq(odeplot(res,scene[i],0..20,thickness=2,color=blue),i=1..3)]));
display(Array([seq(odeplot(res,scenediff[i],0..20,thickness=2,color=red),i=1..3)]));

restart;
Digits:=40;
plot( x*exp(-x),x=0..20,adaptive=true); 
?adaptive
# Quote: The adaptive = geometric option is used for non-parametric plots in the Cartesian coordinate system

The terrible result you get is caused by the default adaptive=geometric.

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