plots:-animate(plot3d,
[[y, x, (1/3)*Pi], x = 0 .. 2*Pi*a, y = 0 .. 5*a, coords = spherical, scaling = constrained,axes=none], a = 0 .. 1);

using spherical coordinates.

d1:=3: d2:=3: d3:=5:    x0:=0: y0:=0: z0:=0:
eq:=[x^2+y^2+z^2-d1^2, (x-3)^2+y^2+z^2-d2^2, x^2+y^2+(z-4)^2-d3^2]:
sph:=[x=x0+r*sin(v)*cos(u), y=y0+r*sin(v)*sin(u), z=z0+r*cos(v)]:
Eq:=simplify(eval(eq,sph)):
R:=min(seq(max(solve(e,r)),e=Eq)) assuming real:
plot3d(rhs~(eval(sph,r=R)), u=0..2*Pi, v=0..Pi, color=green);

int( R^3*sin(v)/3, u=0..2*Pi, v=0..Pi, numeric, epsilon=1e-5); #volume
                          23.55847110
# The integral can be computed symbolically but the result is ugly.

Edit: the code is now a bit more general.

Note first that you are using incorrectly the 2D input. You have a space after assume which is interpreted as multiplication; so all the assume stuff is ignored (actually it is anyway useless in your case).

Maple is not smart enough to obtain such bounds. It does not know Grönwall's inequality (see wiki).

Edit. You should also avoid gamma which is a Maple constant (unless this is what you want).

The series F diverges for any x if n>3.

F:=sum(GAMMA((1/2)*n+1/2)*GAMMA((1/2)*n-1+i)*GAMMA(i+1)^(-i)*((-n+1-2*i)*GAMMA((1/2)*n+1/2)/((x^2-1)*GAMMA((1/2)*n+1/2+i)*(n-1)))^(-i)/(GAMMA((1/2)*n+1/2+i)*GAMMA((1/2)*n-1)), i = 1 .. infinity):
a:=op(1,F):
r:=eval(a,i=i+1)/a:
r1:=simplify(r)/(1-x^2) assuming i::posint:
limit(r1, i=infinity) assuming n>3;

      infinity

For n=3, F = (1-x)/(1+x), 0 <= x < 1.
 

Your procedure newprocedure does not return both matices H and psi; it returns H only.
So, change
return H:
psi:
to
return H, psi:
(actually return is superfluous).

You should also remove interface(rtablesize...) for several reasons, e.g. it is in a loop, you have Dimensions instead of Dimension.

Why don't you save in a text file? E.g.
save newprocedure, "file.mpl";

eq:=Y = (-2*k^3+6*k^2+sqrt(k^8-12*k^7+64*k^6-198*k^5+448*k^4-636*k^3+369*k^2)-7*k-15)/((k^3-3*k^2+5*k-15)*(1+k)):
isolate(eq, indets(eq,sqrt)[])^2:
factor((rhs-lhs)(%)):
select(has,%, Y):
collect(%,Y,factor) = 0;

The equivalence holds if you don't care the branch of the sqrt.

Your code does not work. Try this:

Pairs:=proc(p::prime)
local b, a:=2, t:=0, R:=table():
while (a<p) do
  b := 1/a mod p;
  if isprime(b) and not member([b,a],R)  then t:=t+1; R[t]:=[a,b] fi;
  a:=nextprime(a);
od:
entries(R,nolist)
end:

Pairs(101);
             [7, 29], [37, 71], [43, 47], [53, 61]

Use a simple (linear) change of variables.

J:=2*Int((1/16)*(((1/2)*sin(theta2)*((cos(theta2)+1)*(1+cos(theta1))*cos(phi2*s)+sin(theta2)*sin(theta1))*sin(theta1)/((cos(theta2)+1)*(1+cos(theta1))*(sin(theta1)*sin(theta2)*cos(phi2*s)+cos(theta2)*cos(theta1)+1))))*sin(theta1)*sin(theta2)/Pi^2 * phi2, 
[s = 0 .. 1, phi2 = 0 .. Pi, theta2 = 0 .. Pi, theta1 = 0 .. Pi], method = _CubaCuhre, epsilon = 1e-5):  
evalf(J);
                       0.125000567251692

You must read about how GF(p^k) is represented (e.g. https://en.wikipedia.org/wiki/Finite_field).
If P is an irreducible  polynomial of degree k then a representation for GF(p^k)  is Zp[X] / (P) (i.e. a quotient ring which is actually a field).
So, an irreducible polynomial of degree 4 is needed in your case. For example

G16 := GF(2, 4, 1+x+x^4);

The elements of G16 will depend on this P (but obviously any two G16's will be isomorphic).

You are missing a multiplication sign (* or space) after V[0] and also after a,b,c.
You will need numerical values at least for a,b,c; otherwise the denominator will have to be factored simbolically e.g. with PolynomialTools:-Split and the final result will be huge and probably useless.

You have used square brackets [...]  in eqn1. Replace them with (...)  because [...]  are reserved for lists.

P.S. Spor la treabă!

The best thing is to eliminate yourself the extra vectors due to inherent roundoff errors.
(all vectors after an almost-null one are not reliable).

GS_sel := proc(GS,eps:=1e-7)
  local zer:=false, s;
  s:=proc(v) if zer then NULL elif LinearAlgebra:-Norm(v)>eps then v else zer:=true; NULL fi end;
  map(s,GS)
end:

GS_sel(GramSchmidt([b, c, d, e, a]));