You forgot a multiplication sign * after omega__n
expand(EXPR) works without any assumption.

Just use the Constraints to eliminate variables.

Download Example-vv.mw

# The bug looks severe to me. In Maple 2018 it is OK. 2019?
restart;
fname:= "atest.txt":
for i from 95 to 120 do
  a := 10^i + 13;   aok:=a;
  save a, fname;    read(fname);
  lprint('i'=i, filesize = FileTools:-Size(fname), OK=evalb(a=aok))
od:
# Correct only for i <= 96  (!?)

i = 95, filesize = 107, OK = true
i = 96, filesize = 108, OK = true
i = 97, filesize = 93, OK = false
i = 98, filesize = 94, OK = false
i = 99, filesize = 95, OK = false
i = 100, filesize = 96, OK = false
i = 101, filesize = 97, OK = false
i = 102, filesize = 98, OK = false
i = 103, filesize = 99, OK = false
i = 104, filesize = 100, OK = false
i = 105, filesize = 101, OK = false
i = 106, filesize = 102, OK = false
i = 107, filesize = 103, OK = false
i = 108, filesize = 104, OK = false
i = 109, filesize = 105, OK = false
i = 110, filesize = 106, OK = false
i = 111, filesize = 107, OK = false
i = 112, filesize = 108, OK = false
i = 113, filesize = 109, OK = false
i = 114, filesize = 110, OK = false
i = 115, filesize = 111, OK = false
i = 116, filesize = 93, OK = false
i = 117, filesize = 94, OK = false
i = 118, filesize = 95, OK = false
i = 119, filesize = 96, OK = false
i = 120, filesize = 97, OK = false

Edit. Same for floats! Try:  Digits:=200 and a:=10^i + 13.0

P:=proc(p,x,y)
  local pp:=expand(p), ld:=ldegree(pp,{x,y}), a:=floor(ld/2), b:=ld-a;
  coeff(coeff(pp,x,a),y,b)*x^a*y^b
end proc;

P(13*x^2*y^2 + x*y^2 + 2*y*x^2,x,y) # x*y^2
P(100*x^2*y^2 + 35*y*x + 45*x,x,y)   # 0

The double integral diverges.
I have supposed that Dgamma = D(gamma), Dphi = D(phi).
(Use gamma1 instead of gamma = used by Maple.) 
 

g(a) = piecewise(a < 0, 1/4*2^(1/2)*(a^2)^(1/4)*exp(-1/2*a^2)*BesselK(1/4,1/2*a^2),
-1/2*Pi^(1/2)*2^(1/4)*exp(-1/2*a^2)*(2^(1/2)*CylinderD(3/2,-2^(1/2)*a)+
2*a*CylinderD(1/2,-2^(1/2)*a)));

Note that actually the second branch is valid for a<0 too!

Windows:

src:=cat(kernelopts(mapledir), "\\samples"):
dest:="D:\\tmp\\samples\\":
system( cat("xcopy \"", src, "\" ", dest, " /E") );

No, if   f := y z x + exp(x) cos(y) + g(y, z)
your ff is not correct.

IntWithConst:=proc(f::algebraic, x::name, C::name:=_F)
  local u:=indets(f, And(name, Not(constant))) minus {x};
  int(f, x) + C(u[])
end:

A := exp(x)*cos(y) + y*z:

f:=IntWithConst(A, x);
        f := y*z*x + exp(x)*cos(y) + _F(y, z)

diff(f,y);
       z*x - exp(x)*sin(y) + diff(_F(y, z), y)

The principle is simple: if diff(f(x,y), x) = 0 then f(x,y) is constant with respect to x, which means that f(x,y) depends only on y, i.e. f(x,y) = g(y).

If you want to see this in Maple:

pdsolve( diff(f(x,y), x) = 0, f(x,y) );

         f(x, y) = _F1(y)

Actually, the fact is true only locally: in a non-convex domain, f(x,y) could depend on x (!). Maths is difficult...

So, you want a partial fraction form. Unfortunately Maple does it only for rational expressions and a substitution is needed.

S:=sqrt(N__s):
subs(x=S, convert(subs(S=x,S^2=x^2, ex1), parfrac));

A simpler method is to use the view option

plot3d(4*x^2+9*y^2, x=-5..5, y=-5..5, view=0..100);

The loci are two circles centered at B and C.

(The geometry package cannot be used generally for such problems; they must be solved directly.)

Download geom-loci-ellipse-circles.mw

You have some inaccurate formulations, starting with the definition of a defective eigenvalue (see the cited wiki article).
Anyway, the defective eigenvectors can be obtained directly from the Jordan form of the matrix (they are the columns of Q):

A :=Matrix([[1,-2],[2,5]]):
J,Q:=LinearAlgebra:-JordanForm(A, output=['J','Q']);

v1:=Q[..,1];
v2:=Q[..,2];

(A-3).v1,  (A-3).v2,  (A-3)^2 . v2;

restart;
HB1 := HeunB(5/2, -5^(1/4)*(2*M - R)/sqrt((2*M - R)*(4*M - R)), (17*sqrt(5))/10, (3*5^(3/4)*(2*M - R))/(2*sqrt((2*M - R)*(4*M - R))), 5^(1/4)*sqrt((2*M - R)*(4*M - R))*r^2/(2*(2*M - R)*R^2)):
A:=op(5, HB1):
answer := series(eval(series(eval(HB1, A=z), z), z=A), r);