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. 

What you wrote is far from being the roots of the equation 2*x - tan(x) = 0.
This equation has infinitely many roots (just plot f to see this).

x=1/2*(tan(x))  is an equivalent form of the initial equation.

x=arctan(2x)  is also equivalent but only for -Pi/2 < x < Pi/2;  it has only three roots.

The use of unapply is the best choice.
The first construct does not work because the body of a procedure cannot be altered (without special tools).
A workaround, if you insist in using the -> construct:
seq( (a -> (x -> a*x))(a), a=1..3);

The option view=[...] works as documented. The fact that view=ymin..ymax   also works should be considered as a bonus (and not blamed). Actually for plot3d this is mentioned in the help.

Of course a more organized help system would be nice. In the past, searching the help system was easier and more relevant. It often happens to be more userful a Google search!

- Compute Rank(J) to see whether J is rank defficient.

- If you want to generate matrices with a given rank, multiply a diagonal matrix by invertible matrices. E.g.

with(LinearAlgebra):
A,B:='RandomMatrix(6,6)'$2;
R:=DiagonalMatrix([1$4],6);  #rank=4
J:=A.R.B;
Rank(%);


[A or B could be singular (improbable); then the rank could be <4].

If beta is an integer, the system is polynomial.
But just look at the solution e.g. for beta=4:

solve(subs(beta=4,[f1,f2]),[x,y]);

With Maple there are better and easier to code methods.

restart;
Digits:=15:
ide:=x^2*(diff(y(x), x, x))+50*x*(diff(y(x), x))-35*y(x) - (1-exp(x+1))/(x+1)-(x^2+50*x-35)*exp(x)-int(exp(x*t)*y(t),t=0..1):
n:=10: nx:=16:
p:=x->add( a[k]*x^k,k=0..n):        #  edited
F:=unapply( eval(ide, [y=p]),x):     # (1)  edited
d:=evalf(add(F(k/nx)^2,k=1..nx)):
s:=Optimization:-NLPSolve(d,iterationlimit=20000):
yy:=eval(p(x),s[2]);   # approx solution

 numapprox[infnorm](eval(F(x),s[2]),x=0.1..1); #absolute error for ide
     0.676914366650072e-3

plot(yy,x=0..1);

(1)  Preben noticed an error here; it was essentially   F:=unapply( eval(ide, [y(x)=p(x),y(t)=p(x)]),x):

He also noticed that an exact solution of the ide is exp(x). Let's compare:

numapprox[infnorm](yy-exp(x),x=0..1);
       0.302777540774013e-3