@Cryme 

The procedure uses the general solution obtained by solve (actually SemiAlgebraic).
E.g.
if the solution is  a<x<b, c<y<d then ==>  x = (a+b)/2, y = (c+d)/2
if the solution is  a<x, c<y<x then ==>  x = a+1, y = (c+a+1)/2
etc.

@jnjn0291 

Yes, provided that the eigenvalues are real. Otherwise take the real (or imaginary or both) parts. E.g.

plot3d([seq(Re(R(i)),i=7..8)], x=0..2*Pi,y=0..2*Pi);

@Markiyan Hirnyk 
Just use:
printlevel:=40;  # typesetting=standard

@Carl Love 

The user has m(R+x) (missing `*`). The main problem is using c[1],...,c[4]  and c  followed by  c := ...

@tsunamiBTP 
Yes.
If you just want sinc ant not the DynamicSystems package then define e.g.

sinc := x -> piecewise(x=0, 1, sin(x)/x);

@_Maxim_ 

I totally agree that an epsilon parameter would be nice, but not essential. Note that the epsilon in evalf/Int was also absent in the past. The main problem is the lack of methods.

BTW, the eulermac implementation is far from perfect, see here.

@_Maxim_ 

I don't think that such options are really necessary. Why use an epsilon when Digits should be enough? evalf/Sum usually increases internally Digits.

In your example Maple probably does not have a suitable convergence accelerator.
[The acceleration is usually more efficient for alternating series].

In such cases I think that a combination of symbolic+numeric methods should be used; of course here the symbolic one is sufficient, but take e.g.
evalf(Sum((1-10^(-3))^k/(k+ln(k)), k = 1 .. infinity));

The number of beginning terms is also superfluous because one may use
add(f(n),n=1..N) + Sum(f(n),n=N+1..infinity).

@jacobBN 

With uses in a procedure, the package is not loaded in memory; instead, the necessary functions are invoked by their long names.
with(packagename) can be used only at top-level (i.e. not inside a procedure).

@Markiyan Hirnyk 

Have you read the question or the answer?

@Robert Israel 

Unfortunately there are much simpler expressions for which testeq fails, e.g.

sin((3/7)*Pi)-sin((1/7)*Pi)+sin((2/7)*Pi)-sqrt(7)/2;

@_Maxim_ 

Unfortunately; probably branch problem due to the complex arguments of the elliptics.
MultiSeries seems to be still experimental, so maybe it will be corrected.
 

@_Maxim_ 

From the very beginning:

S1 := int(s(tau),tau=a..t):
E1 := r(t) - diff(r(t),t)*~(S1/s(t)):
MultiSeries:-asympt~(E1,t);

Fine analysis!

I just want to add that the asymptotic behavior of the involute (t-->oo) can be obtained directly using
MultiSeries:-asympt
(for t-->-oo  one has to change t to -t). Of course, MultiSeries:-series(...,t=infinity) is also possible.

@Markiyan Hirnyk 

What should be in your opinion the Maple's answer for your test problem

pdsolve([diff(u(t, x), t, t) = diff(u(t, x), x, x), u(t, 0) = 0, u(t, Pi) = 0]);

Same question in an ideal CAS.