@mehdi jafari 

R gives the parametrization of the surface in spherical coordinates, R = R(u,v). Recall that the Jacobian is r^2*sin(v).
So, the volume equals
int( r^2*sin(v),  r=0..R, u=0..2*Pi, v=0..Pi);
(or equivalently, the double integral above).
 

@brian bovril 

It seems that you have installed a maple init file which does something wrong.
The code works in Maple 2015, 2017 and 2018.

@Carl Love 

Thank you Carl. I think that the situation is clear now.

@Carl Love 

Why are you mentioning only undefined symbols? E.g. Beta is defined and Beta('''f''', A)(x)  acts the same.
The problem is that evalapply is builtin and its action on most functions is not documented. Probably it matters whether the procedure is associated to an operator (such as `[]` and `+`); that's why I have mentioned `&+`.
Note that it's possible to define `evalapply/Beta`  but `evalapply/[]` is ignored.

@Carl Love 

You say "this behavior is a natural and logical result of the overall design".
However, the procedures `[]`, `+`   behave differently than  `&+`, `or`  wrt evaluation of  '''f'''
 

['''f''', A](x);
`[]`('''f''', A)(x);
`+` ('''f''', A)(x);
`&+`('''f''', A)(x);
`or`('''f''', A)(x);
    [('f')(x), A(x)]
    [('f')(x), A(x)]
    ('f')(x)+A(x)
    (`&+`(''f'', A))(x)
    f(x) or A(x)

@Christian Wolinski 

I don't quite understand this type of comments. It would be neater to see the command line followed by its result (given by Maple) and your expected result (+ maybe some explanations). Otherwise it's just some obfuscated code, and the reader is supposed to decrypt it.

On the other side, do you want to program something in Maple and the evaluation of  [...](...)  makes this impossible? In this case please present an example.

@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.