@mehdibaghaee 

You have two elementary mistakes:

1. You define xi := 2*x/a - 1, eta := ...
After that, how could you use xi, eta as new variables in Change?

2. (A more subtle one) You define the procedure U0 in terms of xi, eta (already defined).
This is not correct, you must use unapply in such cases.

@Markiyan Hirnyk 

And what was your "base" recommending  MatrixPolynomialAlgebra:-SmithForm instead of LinearAlgebra:-SmithForm
when they are the same?

@Markiyan Hirnyk 

Have you tried the example?
Please note that exactly the same text appears in the help page of LinearAlgebra[SmithForm].

P.S. You are free to use it for complex coefficients. I will not.

BTW, MatrixPolynomialAlgebra:-SmithForm  calls  LinearAlgebra:-SmithForm
so everything should be clear.

@Markiyan Hirnyk 

It does not work for complex coefficients. Try
A := Matrix(3, 3, [[(27*I)*x+27, 99*x-99*I, 92*x^2-92-(184*I)*x], [(8*I)*x+8, 29*x-29*I, -31*x^2+31+(62*I)*x], [(69*I)*x+69, 44*x-44*I, 67*x^2-67-(134*I)*x]]);

The result should be:

@acer 

Yes, the solution is shorter and cleaner.
My excuse for the longer one is that I prefer to use evala for algebraic numbers/expressions (RootOfs of polynomials) whenever possible.

@brian bovril 

Assume N is odd.
Let A[k] be the vertices on the unit circle of a regular polygon, k=0..N-1.
The line L[k] is tangent to the unit circle at A[k].
The furthest point of intersection is at the distance 1/sin(Pi/2/N) from O
(e.g. intersecting L[0] with L[(N-1)/2]).
The coefficient 1.25 was used to have all the intersections inside the circle with radius R.
a[k], b[k] are the intersections of L[k] with the circle of radius R.
d appears easily if we draw a picture and use some trigonometry.

I obtained the result using solve with series.
Actually also some manipulations by hand because series (and asympt) has some bugs. So, other terms are possible in the expansion:
s =

Edit.

It can be shown that asymptotically:
s = 1/2 * (3*Pi)^(2/3) * r^(1/3)   ≈  2.23 * r^(1/3)

@Preben Alsholm 

It also works for Digits=50.
But for Digits=30 it starts increasing and decreasing Digits and never ends.

Edit. The value for Digits should not be a problem because theoretically Digits will be increased by fsolve if needed.

I never said that Maple has problems only for essential singularities.

To compute a limit could be a very complex task; sometimes almost an "art".
When the expression involves an essential singularity, the result given by Maple has big chances to be wrong.

Unfortunately the help file does not say a word about it, but I am sure that the designers (or at least some of them) are aware of this problem. In these cases, there should appear many "FAIL" answers but this would not be a strong point in marketing.

Here are some real examples (i.e. without a complex direction):

limit((1/2+sin(x))^(1/x), x=infinity);
            1

limit((1+sin(x))^(1/x), x=infinity);
           1

limit( (1+x*cos(x))/(2+x*cos(x)), x=infinity);
           1

@Markiyan Hirnyk 
It seems that you don't understand.
The existence of this f shows that for a CAS it will be very hard to decide whether e.g. the following simple function
g: = sin(exp(-1/z^2+z))*exp(z)-exp(-z)/z^2;
has an essential singularity at z=0 or not.

@Axel Vogt 

The construction of such an f is very technical. It can be done using a deep theorem in complex analysis:
Rosay J.-P., Rudin W. - Arakelian's Approximation Theorem. Amer. Math. Monthly, 96 (1989), 432-434
(the construction itself is not in the article).

@nm 
For a code snippet click the < > icon (near the A icon).

A computation method used by Maple looks very concrete to me.
But the question was not directed exclusively to you.