@ebrahimina 

I don't understand the relevance of your computations for non-square matrices.
If you have two square matrices A,B then the eigenvectors are of course OK.

restart;
with(LinearAlgebra):
m:=4;n:=4;
A:=evalf(RandomMatrix(m,n));
B:=evalf(RandomMatrix(m,n));
Bt := HermitianTranspose(B);
InvBtB := MatrixInverse(Bt . B);
PseudoInvB := InvBtB . Bt;
InvBA := PseudoInvB . A;
ek,vk := Eigenvectors(InvBA)
vkt:=vk^+;

(A-ek[1]*B).vkt[1]^+;  # check ok,  the result is almost 0.

Changing now  m:=4;n:=3;  the result is far from 0.

@Rouben Rostamian  

Nice, vote up!
It would be even nicer to not ignore the overlap of the two cylinders. It seems that Mathematica can do it automatically as a triple integral over their union.
In Maple it is not too difficult. For example, let's compute the mass of the intersection. We take r=1 here.

S:=solve({x^2+z^2<1,  y^2+z^2<1,  x>0},[x,y,z]): # the interesction
remove(hastype,S,`=`);
   [[x < 1, 0 < x, y < -x, -1 < y, z < (-y^2+1)^(1/2), -(-y^2+1)^(1/2) < z], [x < 1, 0 < x, y < x, -x < y, z < (-x^2+1)^(1/2), -(-x^2+1)^(1/2) < z], [x < 1, 0< x, x < y, y < 1, z < (-y^2+1)^(1/2), -(-y^2+1)^(1/2) < z]]

z=-sqrt(-y^2+1)..sqrt(-y^2+1), y=-1..-x, x=0..1: M1:=int(1,%);
M1 := 2/3
z=-sqrt(-x^2+1)..sqrt(-x^2+1), y=-x..x, x=0..1: M2:=int(1,%);
M2 := 4/3
z=-sqrt(-y^2+1)..sqrt(-y^2+1), y=x..1, x=0..1: M3:=int(1,%);
M3 := 2/3
M1+M2+M3; # mass
8/3

So, you have two non-square matrices A,B and want to compute the eigenvalues e_k  for B'A, where B' is the generalized inverse.

What relations do you hope to have between A, B and e_k ? I am skeptical concerning the usefulness of this approach.

@ebrahimina 

If Maple  does not use parallel algorithms for these, it will be very tough for you to implement them.

@ebrahimina 

I don't think that Maple can speed up significantly the computations for this precision.
You will need special tools and hardware. For example there are processors having the floating point arithmetic on 128 bits (==>  34 decimal digits) but it is not widely supported.

P.S. Maybe a better algorithm exists for your problem.

@s144117 

You must be aware that many functions (including hypergeom) are defined (in Maple and other CAS)  beyond the convergence domain (series or integral) using analytic continuation.
A very simple example:
f := hypergeom([1], [], x);
is defined by the series Sum(x^n, n=0..infinity) which converges for |x|<1. But its analytic continuation is f = 1/(1-x) defined for x in C \ {1}.
simplify(eval(f, x=7));
     - 1/6

@s144117 

1. It's not the integral; the series itself diverges (so, the integral is out of the question).

2. It's the ratio test: a series Sum(a[n], n=1..infinity) diverges if Limit(a[n+1]/a[n], n=infinity)  > 1.

restart;
Pn:=proc(n::posint, lambda, mu::Not({0,0.0}), P0, t::name)
local k,P;
P[0]:=P0;
P[1]:=(diff(P[0],t)+lambda*P[0])/mu;
for k to n-1 do  
    P[k+1] := (diff(P[k],t)+(lambda+mu)*P[k]-lambda*P[k-1])/mu
od;
P[n]
end:

Pn(5, 1/4, 1/2, exp(-t), t);
                            - 3 * exp(-t) / 32                               

@Zeineb 

Of course, but only after you define P[0](t), n, lambda, mu.

@sand15 

The OP mentioned "uniform distribution", so probably Prob(X=x) = 1/(k+1).
But anyway he does not seem to be interested in the answer.

I have a technical question, knowing that you are an expert in this kind of things.
1. I have never used thismodule. Trying the example in the help page, it crushed Maple 2018. In Maple 2017 it works. Do you know why? Bug?

2. You have used in this post the construction thismodule['aname']  to access the local variable aname using its global name.
Is this documented?

@student_md 

So, you want to define a sequence of matrices.
It would be better to use as index the position of the matrix in the sequence. I.e. C(k) is the k-th matrix (having the order 2^(k-1) in your case.

You need to define clearly H(k), c(k). With your notation for c_a = c_(m/2) it is not clear what is e.g.  c_3  in the sequence c_0, c_1,...,c_(m/2-1). That's because with your notation, c(m) is defined only for m = "a power of 2".

In other words, use consistent math notations. It will be easy then to obtain the maple code.

@Kitonum 

Of course, but the workaround often works when the locus is not known.
And after all your example is a (degenerate) ellipse.

@Kitonum 

But once we are aware of the (inherent!) problem, workarounds are possible. E.g.
plots:-implicitplot(f(x,y)=d+1e-5, x=0..6, y=0..6, gridrefine=5);

@student_md 

The main advantage is that you can edit it (or even create it with a text editor).

An .m file (internal format) loads faster, but your file is small and does not matter.