To inspect a large matrix without adjusting rtablesize and zooming, you can use the Browse Matrix tool.
Just double-click the placeholder matrix e.g.

A window with scroll bars will open and you will see all the elements.

You should use ApproximateInt only for didactic purposes (or fine tune its many optional parameters).

For example,

Student[Calculus1]:-ApproximateInt(sin(x),x=0..Pi);evalf(%);

    2.008248408

(the exact value being 2).

So, the reliable answers are given by evalf(Int(...))

Without epsilon this could remain unevaluated if the default accuracy cannot be obtained.
(this is because there is a removable singularity at 0 for your integrands).

To "prove" with Maple a nontrivial identity f = g, is is usually not enough.
In most cases, one defines h := f - g; and try to simplify to 0 using all sorts of acrobatics with simplify, convert, expand, assuming etc.

maximize(F,location);

works, as tomleslie stated.

Note that simplify(F) gives a better reprezentation.

Note also that your function is not continuous. If you use the almost equivalent definitition

G:=piecewise(t<0,0,t<=1,t^2,t<3,3-t);

then maximize(G,location) gives the same result

    2, {[{t = 1}, 2]}

but t=1 is not actually the maximum point, because the max is not attained.

Unfortunately the default method gives a wrong result. The methods ftoc and ftocms work correctly.

int(1/(x^4+4), x=0..1, method=ftoc);

    1/16*ln(5)+1/8*arctan(2)

Maplesoft should do something if such a simple integral is wrong!

The symmetry exhibited here is due to the real analyticity of F i.e. F(s) is real for a real s and so
conjugate(F(s)) = F(conjugate(s)).

It's not only the complexity. Your huge integrand contains unknown (arbitrary) functions and their derivatives. What kind of result do you expect to obtain?

Maple uses here wrongly some "generic" formulas.

A simpler version is:

Sum(cos(8*Pi*j/n),j=0..n-1);

value(%);
      0

This is true for all positive integer n except  n = 1  or 2  or  4.

You have a simple access to each pixel (using an Array), so you can apply any transform.
If the procedure is compiled, it will be fast.
See an example here

A very simple one:

with(ImageTools):
X:=Read(cat(kernelopts(mapledir),"/data/images/fjords.jpg")):
dims:=rtable_dims(X):
m,n:=rhs(dims[1]),rhs(dims[2]):
f:=(i,j,k)->X[round((i-1)^3/m^2)+1,j,k]:
Y:=Array(dims,f,datatype=float[8]):
Embed([X,Y]);

Just use:

abs~([zeros]);

The root with largest absolute value:

sort([zeros],key=abs)[-1];

You may use polar coordinates like this:

eval(x-a+y-2, [x=a+rho*cos(theta), y=b+rho*sin(theta)]):
limit(%, rho=0);
                             b - 2

But be aware that this is not exactly the same thing, because the previous limit is not uniform in theta
(i.e. one obtains this way only a directional limit, and the actual limit(..., {x=a,y=b}) may not exist
-- not in this case).

It seems that the Psi expressions are difficult to simplify.
But noticing (as in your link) that our integral J(n) is the n-th derivative wrt a of

Int( x^(a)/(x^3 + 1), x=0..infinity);

for a=0, it is easy to compute J(n):

Ja:=int( x^(a)/(x^3 + 1), x=0..infinity):
VJa:=value(eval(diff(Ja,a$n),a=0)):
J:=unapply( 'simplify'(VJa), n):

seq('J'(n)=J(n),n=1..10);

Edit. This also works for the integral

for p>m+1, n>=0, m>-1. Simply replace
Ja:=int( x^(a)*x^m/(x^p + 1), x=0..infinity):

I would not trust it except for extremely simple expressions.

Some examples:

restart;
assume(t>-13);
f:=t^4-(t+1)^4-(t-2)^4;
                     4          4          4
                    t  - (t + 1)  - (t - 2)
f:=simplify(%);
                   4      3       2            
                 -t  + 4 t  - 30 t  + 28 t - 17
coulditbe(f>0);  #?
                              true
coulditbe(f<0);  #?
                              FAIL
is(f>0);
                             false
is(f<0);
                              true
coulditbe( x^7+y^7=z^7) assuming posint;  #?
                              true
coulditbe( x^n+y^n=z^n) assuming posint, n>2;  #why, if previous is true?
                              FAIL

P.S. Of course we must accept a FAIL answer from is and coulditbe. But a wrong true/false is a bug.

select(u->is(u>0),  [solve(%, omega)]) assuming positive;

dsolve(ode)  finds it. It is given implicitely.
You may use dsolve(ode, explicit)  but an inherent RootOf will appear.