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]));

Just add

dsolve({ entries(de, nolist), x__1(0)=0, x__2(0)=10, D(x__1)(0)=0, D(x__2)(0)=0 }) ;

Note that I have included initial conditions for derivatives, otherwise two arbitrary constants appear.

The minimal correction would be to replace both rand(1..6);  with rand(1..6)();
That's because rand(1..6)  is a procedure and must be called to produce a random number in the interval 1..6.

It will work, but this way each time you call it, the two procedures are created and then called.
A more efficient way is to define
dice := rand(1..6), rand(1..6);   # or better  rand(1..6) $2

Now dice is no longer a procedure but acts exactly the same, i.e. 
dice() ;
       6, 4

Try this one

COEFF:=(f, z) -> `if`(has(z,`&*`), coeff(args), coeff(eval(f,`&*`=0),z,_rest));