@Markiyan Hirnyk  Please explain how it could happen that in the post on January 10 2016, all comments have been removed, including mine, this post has turned into a question with your answer in the form my animation?

You wrote in brackets something like this

restart;

x=solve(y, x);

                                    x=()

@vv  Why to check, we can not instead of

eval(  'evalf'( 'add'('add'(a, i=1..N), j=1..N)),  N=1000);

simply write

add(add(evalf(a, i=1..1000), j=1..1000);

The result is the same, but it requires 10 times less time. 

@Majmaj You can't use  labels = ["x values", "y values"] if axis[2] = [color = white], because everything that refers to  y-axis  will be white.

Markiyan Hirnyk  Why do you never specify the source of the problem? This is very similar to the problems of student competitions.

@tomleslie  for the adjustment. I did not know before about this option.

Unfortunately, some old bugs remain. Here is the example of the definite integral with assuming:

A:=int(abs(cos(n*x)), x=0..Pi);

A assuming n=0;

eval(A) assuming n=0;

value(A) assuming n=0;

evalf(A) assuming n=0;

simplify(A) assuming n=0; 

Of course the integral is easy to calculate by the direct substitution:

eval(A, n=0);

I have put  +1 for OP,  but my vote up inexplicably disappeared.

Addition: Now I see that it appeared again.

Good work, vote up. Just one remark. Your code is too complex for this problem.  Roughly the same can be done much more simple means:

A := plottools[polygon]([[sqrt(3)+1, 0, 0], [1, 1, 0], [1, -1, 0]], color = "Violet"):

B := plottools[polygon]([[sqrt(3)+1, 0, 0.01], [1, 1, 0.01], [1, -1, 0.01]], color = "Green"):

C := plots[display](A, B):

F := t->plottools[rotate](t, alpha, [[0, 0, 0], [0, 0, 1]]):

G := t->plottools[rotate](t, alpha, [[0, 0, 0], [0, 1, 0]]):

plots[animate](plots[display], [(G@F)(C)], alpha = 0 .. 2*Pi, scaling = constrained, axes = none, frames = 180, lightmodel = none);

 The idea of the code - the triangle at the same time rotates around the y-axis and z-axis.

Addition. A slight modification of the code clearly shows that the rotations in 3d around different axes do not commute:

A := plottools[polygon]([[sqrt(3)+1, 0, 0], [1, 1, 0], [1, -1, 0]], color = "Violet"):

B := plottools[polygon]([[sqrt(3)+1, 0, 0.01], [1, 1, 0.01], [1, -1, 0.01]], color = "Green"):

C := plots[display](A, B):

F := t->plottools[rotate](t, alpha, [[0, 0, 0], [0, 0, 1]]):

G := t->plottools[rotate](t, alpha, [[0, 0, 0], [0, 1, 0]]): 

Anim1:=plots[animate](plots[display], [(G@F)(C)], alpha = 0 .. 2*Pi, scaling = constrained, axes = none, frames = 180, lightmodel = none):

Anim2:=plots[animate](plots[display], [(F@G)(C)], alpha = 0 .. 2*Pi, scaling = constrained, axes = none, frames = 180, lightmodel = none):

plots[display](<Anim1 | Anim2>);

@kountouretti  For visualization I used my code in the answer below  "Direct calculation". I took  epsilon = 0.1. We find that for all  n>=287  all the points  (n, x(n))  lie within the grey band. To the point of not merged into one line I took  n  with the step = 10 :

restart;

x := n->n^2/(n^2+31*n+228):

x0 := limit(x(n), n = infinity):

N := epsilon->max(1, floor(op(1, solve(abs(x(n)-x0) <= epsilon, n)[-1]))+1):

 N(0.1);

                                                      287

A := plot([1-.1, 1+.1, [seq([n, x(n)], n = 1 .. 500, 10)]], n = -2 .. 500, -0.1 .. 1.25, style = [line, line, point], color = [black, black, red], linestyle = [3, 3], symbol = solidcircle, symbolsize = 7, labels = [n, "x(n)"], tickmarks = [default, [0.9, 1, 1.1]]):

B := plots[inequal]({y > 0.9, y < 1.1}, x = 0 .. 500, y = 0 .. 1.25, optionsfeasible = [color = "LightGrey"]):

plots[display](A, B);

Edited. 

@AdVen1935   See the corrected file. To upload a file click the thick green up-arrow in the menu of the editor.

Histogram1.mw

@AdVen1935  See attached file. I was a little changed procedure code and your example code.

Histogram.mw

@Preben Alsholm  

Thanks for the detailed explanation.

I decided to once again carefully check everything. Your version works in a Standard Interface 32 bit Maple 2015.2 (under Windows 8) when we first start. But if, for example, remove the line of code  Expr := 'Expr'; , and then re-insert the code stops working. In Standard Interface  Maple 16  your code does not work.

In the Classic Interface  Maple 2015.2  your code (and my last version too) works without any limitation.

@Preben Alsholm  Brilliant solution!

Here's another version ( slightly shorter), based on the same idea:

Expr:=a*diff(x(t),t)+b*x(t)+r*diff(x(t),t,t)+a*diff(y(t),t)+b*diff(y(t),t,t)+c*y(t);

S1:=op~([selectremove(s->has(s, diff), [op(Expr)])]);

Expr:='Expr':

gc():

Expr:=add(i, i=S1);

Unfortunately your code (and mine too) does not work in the standard interface Maple 2015.2.  Only in classic one.

@Jenser  Perhaps you have an older version of Maple. Replace the line of code 

``(add(select(s->has(s,diff), L)))+add(select(s->not has(s,diff), L));

by the line

``(`+`(op(select(s->has(s,diff), L))))+`+`(op(select(s->not has(s,diff), L)));