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

@Bendesarts  The equation  eq_liaison  can be simplified by one-lined code:

eq_liaison:=sin(gamma0(t))*cos(beta0(t)) = -(sin(gamma[1](t))*sin(psi[1](t))*sin(theta[1](t))-sin(gamma[1](t))*cos(theta[1](t))*cos(psi[1](t))+cos(gamma[1](t))*sin(psi[1](t))*cos(theta[1](t))+cos(gamma[1](t))*cos(psi[1](t))*sin(theta[1](t)))*cos(beta[1](t));

lhs(eq_liaison)=1/2*(trigsubs(2*combine(rhs(eq_liaison))))[];

I do not understand your desire maximally to automate the task of fairly particular. I believe that the most effective way to work in Maple - is a reasonable combination of manual and automatic operation. The criteria of effectiveness - is the most quick and easy solution to a specific problem.

@Carl Love  Good solution! Unfortunately, I do not know much English and for me often easier to write new code than to understand the capabilities of the existing built-in commands.

@vv  I agree, but I understand the OP's question as finding the exact value  N(epsilon)

@Age  The quotes in  'deriv'(x)  exclude the premature  calculation. Without quotes  x  in  deriv(x)  should be a constant.

See

deriv(x);

Error, (in deriv) invalid input: fdiff expects its 2nd argument, N, to be of type {integer, name, list(integer), list(name), list(name = constant), set(name), set(name = constant), name = constant}, but received t = x