Maple can approximate the integral (you can use evalf), but cannot compute it symbolically because the antiderivative cannot be obtained in termes of Maple functions.
The idea is to use the residue theorem.  In your case:

restart;
f := (a,z) -> z*exp(a*z) / (z^2 + 1)^2:
R:=[singular(f(a,z),z)];
2*Pi*I*add(residue(f(a,z),eq[]), eq=R): simplify(%);

                    R := [{z = -I}, {z = I}]
                         I Pi a sin(a)

(The result is valid for any Jordan curve positively oriented an surrounding the poles).

See wiki.

stp:=L -> [seq([L[i], L[i]+~L[i+1]][], i=1..nops(L)-1), L[-1]]:
L2f:=u -> `if`(u[2]=0,infinity,u[1]/u[2]):
Brocot:= n -> L2f~((stp@@n)([[0,1],[1,0]])): 

Vector([seq(B||n = Brocot(n),n=0..5)]);

(You may want 1/`0` instead of infinity).

BesselJ(3+1/2, x) simplifies to an elementary expression

and catastrophic cancellation occurs for small x.

For BesselJ(3.5, x) a numeric approximation is used and the the loss of significance is taken into account by the algorithm.

You can increase Digits.

The problem is that you use the variable Zeta which is actually a (celebrated) Maple function.
Change its name to ZETA.

Ec := (Ems+I*Eml)*(1+((Ems+I*Eml)/Ef-1)*Zeta*phi/((Ems+I*Eml)/Ef+Zeta))/(1-((Ems+I*Eml)/Ef-1)*phi/((Ems+I*Eml)/Ef+Zeta)):
EC:=eval(Ec,Zeta=ZETA):
simplify(evalc([Re,Im](EC)));

Do not use [...] for grouping in expressions; they are reserved for lists.
Do not use floats in symbolic computations.

restart;
f:=1/(exp((-1+(1-4*epsilon)^(0.5))/(2*epsilon))-exp((-1-(1-4*epsilon)^(0.5))/(2*epsilon)))*(exp((-1+(1-4*epsilon)^(0.5))/(2*epsilon)*t)-exp((-1-(1-4*epsilon)^(0.5))/(2*epsilon)*t))/(exp(1-t)-exp(1-t/(epsilon))):
g:=convert(f, rational);

simplify(limit(g, epsilon=0, right)) assuming t>0;
                               1
simplify(limit(g, epsilon=0, left)) assuming t>0;
                               0
simplify(limit(g, epsilon=0)) assuming t>0;
                           undefined

So, for t>0 your function has distinct side-limits. (For t=0 the function is not defined).

Download puiseux.mw

SymFun:=proc(f,x)
local s,a_,b_, g;
if is(simplify(diff(f,x))=0) then return 'ConstantFunction' fi;
s:=solve(identity(eval(f,x=a_+x)=eval(f,x=a_-x),x),[a_]);
if s<>[] then return 'VerticalSymmetry'=s[] fi;
s:=solve(identity(eval(f,x=a_-x)+eval(f,x=a_+x)=2*b_,x),[a_,b_]);
if s<>[] then return 'CentralSymmetry'=s[] fi;
# hack
g:=eval(f,x=a_+x)-eval(f,x=a_-x);
s:=solve(limit(g,x=infinity),[a_]);
if s<>[] and is(simplify(g,s[1])=0) then return 'VerticalSymmetry'=s[1] fi;
'NoSymmetry'
end:

eliminate({s = sum(n/2^k, k=0..p), n/2^p=1}, p):
simplify(solve(%[2],s))[];
              s = -1+2*n

A point on the unit sphere is determined by the spherical coordinates (th, ph). You can choose any color for this point. E.g.

c := ph+th/16:  # color (HUE)
plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical, color=c);

If your vector field is V(r,th,ph), you can take c := LinearAlgebra:-Norm(V(1,th,ph), 2, conjugate=false);

Digits:=15:
par:=[lambda1 = 0.733e-1,lambda2 = 5.3344,alpha = 4.8492]:
f := (Z^(lambda2/lambda1))^2*ln(Z)^2/((Z-1+alpha)^2*(Z^(lambda2/lambda1)-1+alpha)^3):
J:=IntegrationTools:-Change(Int(f,Z=1..infinity), Z=z^(lambda1/lambda2)) assuming positive:
evalf(eval(J,par));

           1.15034205858649*10^(-7)  

Download modular-cp.mw

restart:
A:= <x1,y1,z1>:
B:= <x2,y2,z2>:
P :=<x3, y3, z3>:
H:=A+t*(B-A):
solve( (P-H)^+ . (B-A), t ):
'H'=eval(H, t=% );