BIN:=module()
export ModuleApply, Fac, Fac1;
local pmax:=trunc(2^((kernelopts(wordsize)-1)/2));

Fac:=proc(n,m,p)
  local r:=1,k;
  for k from n to m do r := irem(r*k,p) od
end:

Fac1:=proc(n::integer,m::integer,p::integer)::integer;
  option autocompile;
  local r::integer:=1, k::integer;
  for k from n to m do r := irem(r*k,p) od
end:

ModuleApply := proc(N::nonnegint, K::nonnegint, p::prime)  # combined with Carl's
local r:=1, n:=N, k:=min(K, N-K), Bin, fac:=`if`(p<pmax,Fac1,Fac);

Bin:=proc(N,K)  # 0 <=K < p, 0 <= N
local k:=min(K,N-K);
if k<0 then return 0 fi;
fac(N-k+1,N,p)/fac(1,k,p) mod p
end;

if k<0 then return 0 fi;
while k > 0 and r > 0 do 
   r:= irem(r*Bin(irem(n,p,'n'), irem(k,p,'k')), p) 
od;
r
end:

end module:

p:=nextprime(8*10^7): n:=p*2-1:  k:=floor(p/2):
CodeTools:-Usage(BIN(n,k,p));
memory used=2.97MiB, alloc change=0 bytes, cpu time=1.89s, real time=1.88s, gc time=0ns
                            80000022

p:=nextprime(9*10^8): n:=p*2-1:  k:=floor(p/2):
CodeTools:-Usage(BIN(n,k,p));
memory used=2.11KiB, alloc change=0 bytes, cpu time=20.78s, real time=20.78s, gc time=0ns
                           900000010

NOTE. The compiler works for p < 3.037*10^9. For p greater than that, BIN uses the non-compiled Fac but then k should be smaller (or at least the reductions via Lucas' theorem) to have reasonable execution time. The magnitude of N is not important.
 

macro(ff=2*a*b-3*a-2*b-c):
f:=(a,b,c)-> `if`(ff<=30 and ff>=15,ff,0):
A:=Array(0..10,1..5,0..1,f, storage=sparse):
indices(A,pairs); # or simply inspect the Array using the matrix browser

(10, 3, 1) = 23, (7, 4, 1) = 26, (5, 4, 0) = 17, (4, 5, 1) = 17, (9, 3, 0) = 21, (8, 3, 1) = 17, (4, 5, 0) = 18, (6, 4, 0) = 22, (5, 4, 1) = 16, (8, 3, 0) = 18, (10, 3, 0) = 24, (5, 5, 0) = 25, (7, 4, 0) = 27, (5, 5, 1) = 24, (9, 3, 1) = 20, (7, 3, 0) = 15, (6, 4, 1) = 21

Vector needs a procedure (with 1 argument) as initializer.
fenq  is such a procedure but  fsenq(2.1, j)   (actually you mean fseqn) is not; it's an expression.
So, use

psf := Vector(5, j->fseqn(2.1, j));

or maybe the more "sophisticated"

psf := Vector(5, curry(fseqn, 2.1));

The roots of integer numbers are not automatically simplified. Not even 4^(1/2). It is by design.
Unfortunately it is hard to find in the help pages "automatic simplification",  even if the Programming Guide contains information about this important aspect.

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]