It is enough to change:

Num_N :=simplify(convert(series(Num, u, N),polynom)):  #
Den_N :=simplify(convert(series(Den, u, N),polynom)):  #

S:= n -> evalf(Sum((1/k^0.1)*sin(1/k), k=1..n)):
S(1000000);

                          7.883752545

For x0 = 0 you have y0 = exp(-1) and the correct series is:

Edit. This is for the equation ln((1+x)*y) + exp(x^2*y^2) = 0
which appears first.
For ln((1+x)*y) + exp(x^2*y^2) = x + cos(x)   see:
https://www.mapleprimes.com/questions/223411-Real-Function-Defined-By-A-Implicitly-Equation

It would be useful to know (for the future) how did you manage to introduce that special Unicode character "∷"  (named "proportion").  Copy & paste from a pdf file?

It is difficult to debug Maple code containing such characters.

You should never unprotect and redefine Pi.
(unprotect is designed for advanced users only).
Just use pi, not Pi.

Groebner:-Basis actually computes the reduced Groebner basis.

It is disappointing.
But at least for floats it is correct:

conic(cf, evalf(eq), [x, y]):
coordinates(center(cf));
       [.5000000000, -.8660254038]

The eigenvectors are associated to eigenvalues, which are not sorted (they appear in a random order).

You have the possibility to sort the eigenvalues and pick the corresponding eigenvectors.

Change the unknown function f:

F(tau) = int( exp(lambda*t)*f(t), t=0..tau ); hence f(t) = F'(t)*exp(-lambda*t), F(0)=0

so, the problem reduces to an ODE in F.

This works:

restart;
with(Statistics):
DD := Distribution(ProbabilityFunction = (k -> piecewise(k > 0, 2^(-k))), 'Support' = 1 .. infinity):
# ProbabilityFunction(DD,t);
Mean(DD), Variance(DD);

        2, 2

aa is not a polynomial; x^(1.25) ?

On my computer fsolve finishes in 25 sec, Analytic in 30 sec (Maple 2017, 64 bit).
(maxsols is superfluous for non-polynomial functions).

Analytic uses special algorithms but works only for analytic (=holomorphic) functions; fsolve uses a Newton-type method.
E.g.  
fsolve(abs(sin(z))-z/10, z=5..20);   # works, but 

RootFinding:-Analytic(abs(sin(z))-z/10, z=5-I..20+I/10);  # does not.

Probably the 2d parser translates

2.59/3.7, 3.42/3.8

    .7000000000, .9000000000

to:

2.59*(1/3.7), 3.42*(1/3.8);

    .7000000001, .8999999999

SolveTools:-SemiAlgebraic({c*x^7-6*x^3+6=0}, parameters=[c]);

piecewise(c < 0, [[x = RootOf(_Z^7*c-6*_Z^3+6, index = real[1])]], c = 0, [[x = 1]], c < RootOf(823543*_Z^3-1492992, 1.188 .. 1.25), [[x = RootOf(_Z^7*c-6*_Z^3+6, index = real[1])], [x = RootOf(_Z^7*c-6*_Z^3+6, index = real[2])], [x = RootOf(_Z^7*c-6*_Z^3+6, index = real[3])]], c = RootOf(823543*_Z^3-1492992, 1.188 .. 1.25), [[x = RootOf(50421*_Z^5*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)^2+62208*_Z^6+74088*_Z^4*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)-117649*_Z^2*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)^2+108864*_Z^3-172872*_Z*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)-254016, index = real[1])], [x = RootOf(50421*_Z^5*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)^2+62208*_Z^6+74088*_Z^4*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)-117649*_Z^2*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)^2+108864*_Z^3-172872*_Z*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)-254016, index = real[2])]], RootOf(823543*_Z^3-1492992, 1.188 .. 1.25) < c, [[x = RootOf(_Z^7*c-6*_Z^3+6, index = real[1])]])

P.S. It could be instructive to solve it by hand; it is not difficult, writing the equation as  c = f(x).

Edit. If you want only positive roots, use:

SolveTools:-SemiAlgebraic({c*x^7-6*x^3+6=0,x>0}, parameters=[c]);

Replace:

plot(x, x=0..1);
with
print(plot(x, x=0..1));

(this is because of how printlevel works, see ?printlevel).