Yes, it seems that

seq(C[i+1]-C[i], i = 1 .. N-1):
does not finish the job.
But,
seq(C[i+1]-i^2, i = 1 .. N-1):
works.
Probably a bug.

The limit does exist in C.
For a real function it is possible to use

RealDomain:-limit(sqrt(x^2 - 3 *x + 2), x=2, left);
        undefined

Or, execute first
with(RealDomain);

However I'd recommend RealDomain only if the user really does not know (yet) about complex numbers & functions.
(You should know that almost all mathematical functions are defined in C.)

FC:= proc(n) 
uses plots,plottools;
local a,b,u;
{seq(seq(a/b, a = 0 .. b), b = 1 .. n)};
display(seq(circle([u,1/2/denom(u)^2],1/2/denom(u)^2), u=%)
        ,scaling = constrained, color = blue, size = [800, 800])
end proc:

FC(6);

See also Apollonian circles

print~([a, b]);  means  [print(a), print(b)];

Now, print(a)  prints a (as a side effect) and returns NULL.
The last output is [NULL,NULL]. If you want to eliminate it, use print~([a, b]):  (colon instead of semicolon).

Digits:=150:
A:=Matrix(eval(M,P=0),shape=symmetric):
B:=Matrix(-diff(M,P), shape=symmetric):
E:=LinearAlgebra:-Eigenvalues(A,B):
E1:=Re(E);
E2:=Im(E); # =0
Esorted:=sort(E1);
evalf[15](Esorted[1..15]), evalf[15](Esorted[-15..-1]);

Yes, shape=triangular[upper, unit]  seems to have a serious problem.

The best (partial) workaround is probably

U:= Matrix(Matrix(4,4, shape=triangular[upper, unit], symbol =u), shape=triangular[upper]);

(partial, because the diagonal elements are now assignable.)

s:=solve((hv*hw+lv*lw)/(lv+hv)*(lv+hv) = dw*(lv+hv), hv):
sign(s)*numer(s)/sort(expand((sign(s)*denom(s))),order=plex(hw,dw));
#   Do you think it's worth the effort?

solve({seq(a*v + b*u - w)}, {a,b});

(If w is not in the span of u,v, the result is NULL).

Replace n:=10^10 by  n:=10.^10

Actually, the problem is about inequalities in solve, in your case

F := 1-(1/2)*t*exp(-(1/2)*t)-exp(-(1/2)*t):
S := solve({F=u, t >=0}, t) assuming u >= 0, u <= 1;

In general, solve is able to manage inequalities only for polynomials (of several variables).
In most other cases, the inequalities are simply ignored. So, here, S = S1, where
S1 := solve({F=u}, t) assuming u >= 0, u <= 1;

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