The "end do" in the loop should clearly be at the end.
So with that corrected your code runs:
 

restart;
with(plots):
with(plottools):
K := 48;
R1 := 30; R2 := 20; OA := R1+R2;
Vodilo := PLOT3D(cuboid([-3, -3, 0], [OA+4, 4, 2]));
Koleso1 := PLOT3D(cylinder([0, 0, 5], R1, 18), cylinder([0, 0, -2], 1, 16));
Koleso2 := PLOT3D(cylinder([0, 0, 7], R2, 6), cylinder([0, 0, -2], 1, 16));
Amplituda := .3:
###
for i from 0 to K do   
  phi := sin(6.28*i/K)*Amplituda;
  phi1 := phi*(R1+R2)/R2;
  P1 := rotate(Koleso2, 0, 0, -phi1);
  P2 := rotate(Vodilo, 0, 0, -phi);
  x := OA*cos(phi); 
  y := OA*sin(phi);
  P[i] := display(Koleso1, translate(P1, x, y, 0), P2)
end do: # colon
###

display(seq(P[i], i = 0 .. K), orientation = [55, 101], insequence = true, scaling = constrained, axes = normal);

This unassigns all the numbered members:
 

restart;
with(Statistics):
X := a -> RandomVariable(Exponential(a)):
##
for n from 1 to 100 do
  S := X(n);
  # here are some operations to do
end do:
##
A:={anames(alluser)}:
IA:=indets(A,suffixed(_ProbabilityDistribution,integer)):
##
eval(_ProbabilityDistribution5);
exports(_ProbabilityDistribution5);
##
unassign(op(IA));
eval(_ProbabilityDistribution5);

But you might want to unassign all but the last one:
 

N:=100:
for n from 1 to N do
  S := X(n);
  # here are some operations to do
  if n = N-1 then 
     A:={anames(alluser)};
     IA:=indets(A,suffixed(_ProbabilityDistribution,integer));
     unassign(op(IA));
  end if
end do:
anames(alluser);

Your integral equation has infinitely many solutions because the corresponding homogeneous equation has a one parameter family of solutions.
A particular solution can be found as vv did.

You will find my reasoning in the attached worksheet:

Download IntegralEquations_delay_Z.mw

I have corrected some syntactical errors. But dsolve cannot handle this analytically.

Certainly numerically, though:
 

restart;
ode:=m*diff(x(t),t,t)= - m * A * sin(2*Pi*f*t) - piecewise( abs(x(t)) < x__max, 0, abs(x(t)) >= x__max, -k*(abs(x(t))-x__max)*signum(x(t))) - diff(x(t), t);
dsolve(ode); #NULL
nms:=indets(ode,And(name,Not(constant))) minus {x,t};
res:=dsolve({ode,x(0)=x0,D(x)(0)=x1},numeric,parameters=[nms[],x0,x1]);
res(parameters=[A=1,m=1,f=1,x__max=2,k=2,x0=0,x1=1]); # Example
plots:-odeplot(res,[t,x(t)],0..20);

Try this:

restart;
Rok := Vector([2013, 2014, 2015, 2016, 2017, 2018], datatype = float);

TrzbyCelkemEmco := Vector([1028155, 1134120, 1004758, 929584, 995716, 1152042], datatype = float);
########
KubickaTrzby := Statistics:-PolynomialFit(3, Rok, TrzbyCelkemEmco, x);
p1:=plot(KubickaTrzby,x=2013..2018):
p2:=plot(Rok,TrzbyCelkemEmco,color=blue):
plots:-display(p1,p2); # Not impressive
########
pol:=a*x^3+b*x^2+c*x+d;
resid:=eval~(pol,x=~Rok)-TrzbyCelkemEmco; # Residuals
A,B:=LinearAlgebra:-GenerateMatrix(convert(resid,list),[a,b,c,d]);
op(1,A); # Dimension of A: 6x4
LinearAlgebra:-Rank(A); # 2 only

The rank could in principle have been min(6,4) = 4. That would have meant full rank.

To get an exact answer to this question is certainly not easy to say the least.

But the beginnings of an attempt is this:
 

restart;
### Your right hand sides:
fS := -D2*x1*x2-D1*x1-x1*x3+S; 
fV := D2*x1*x2-D4*x2+x1*x3; 
fC := -D6*x3*x4-D5*x3+x2; 
fR := D6*x3*x4-x4;
eqs:=solve({fS,fV,fC,fR},{x1,x2,x3,x4}); # Notice RootOfs
L:=map(subs,[eqs],[x1,x2,x3,x4]); # The equilibrium points
### The Jacobian matrix as a function of x1,x2,x3,x4:
J:=unapply(VectorCalculus:-Jacobian([fS,fV,fC,fR],[x1,x2,x3,x4]),x1,x2,x3,x4);
J(x1,x2,x3,x4); # View
J(op(L[1])); # Point number 1
EV1:=LinearAlgebra:-Eigenvalues(%);
### The first two are negative since D1>0.
### But the two other depend on the actual size of your parameters as seen here:
### For these values the system is asymptotically stable:
eval(EV1[4],{D1=0.1,D2=1,D4=0.1,D5=.1,S=.0001});
eval(EV1[3],{D1=0.1,D2=1,D4=0.1,D5=.1,S=.0001});
### For these values the system is unstable:
eval(EV1[4],{D1=0.1,D2=1,D4=0.1,D5=.1,S=1});
eval(EV1[3],{D1=0.1,D2=1,D4=0.1,D5=.1,S=1});

This just illustrated the problem for point 1.

Without specifying Q concretely you get:
 

restart;
ode:=diff(v(x),x,x)+diff(v(x),x)*cot(x)-v(x)*(cot(x)^2+nu)=a^2*Q(x)/D1;
dsolve(ode);

Notice that there is an unevaluated integral involved simply because Q is not known.

That NULL is returned means that solve couldn't find a solution or that none exists.

In fact the latter is the case here with your parameters so far just names.
Your system is linear in the a-variables, the coefficient matrix A is 9x9 and its rank is only 8.
Thus for general right hand side B in A.x = B there is no solution for x = <a[0],..., a[8] >.
 

restart;

omega := v/h;
t := sum(a[j]*x^j, j = 0 .. 6)+a[7]*cos(omega*x)+a[8]*sin(omega*x);
r1 := diff(t, x$2);
r2 := diff(t, x$4);
c1 := eval(t, x = q+2*h) = y[n+2];
c2 := eval(r1, x = q) = f[n];
c3 := eval(r1, x = q+h) = f[n+1];
c4 := eval(r1, x = q+2*h) = f[n+2];
c5 := eval(r1, x = q+3*h) = f[n+3];
c6 := eval(r2, x = q) = g[n];
c7 := eval(r2, x = q+h) = g[n+1];
c8 := eval(r2, x = q+2*h) = g[n+2];
c9 := eval(r2, x = q+3*h) = g[n+3];
b1 := seq(a[i], i = 0 .. 8);
#####
A,B:=LinearAlgebra:-GenerateMatrix([c1, c2, c3, c4, c5, c6, c7, c8, c9],[b1]);
LinearAlgebra:-Dimension(A); #9x9
LinearAlgebra:-Rank(A); # 8

Use fsolve instead. Your function is strictly increasing on the real axis, so the equation f(x) = 5 has at most one real solution.
It is odd and tends to infinity as x -> infinity. Thus there is exactly one real solution.

fsolve(f(x)=5,x);
############## Workaround ############

restart;
f := x -> x*sqrt(4*x^2 + 1) + arcsinh(2*x);
convert(f(x),ln);
solve(%=5,x);
allvalues(%); # Only one
evalf(%);
#### With 5.0 instead of 5:
restart;
f := x -> x*sqrt(4*x^2 + 1) + arcsinh(2*x);
convert(f(x),ln);
solve(%=5. ,x);

You obviously by pi mean the mathematical constant, which in Maple is Pi.

Ideally, you would just do:
 

restart;
ode := diff(y(x), x, x) + (n*Pi)^2*y(x) = A^3*sin(n*Pi*x)^3;

dsol1 := dsolve({ode, y(0) = 0, y(1) = 0}) assuming n::integer ;

It results in a division by zero error.

Are you sure that there is a twice differentiable solution?
Try letting n and A both be 1 for a concrete example:
 

dsolve({eval(ode,{n=1,A=1}), y(0) = 0, y(1) = 0});

No solution is found.
Worse:
 

ode1:=eval(ode,{n=1,A=1});
sol:=dsolve({ode1, y(0) = 0, D(y)(0) = y1});
eq:=eval(rhs(sol),x=1);

Since eq = 3/(8*Pi) is independent of y1, no value of y'(0) = y1 will make y(1) = 0.

Thus there is no solution to your problem (here only shown for n=1 and A=1).

A general right hand side for n=1:
 

restart;
ode := diff(y(x), x, x) + Pi^2*y(x) = f(x);
sol:=dsolve({ode,y(0)=0,D(y)(0)=y1});
eq1:=eval(sol,x=1);
## Examples:
value(eval(eq1,f=sin)); # not zero
eval(eq1,f=sin^2);
value(%); 
evalf(%);

You need f to satisfy: int(sin(Pi*x)*f(x), x = 0 .. 1) = 0.

A simple example of a left hand side for which there is a solution is f(x) = x - 1/2, because eq1 is satisfied:
 

ode1:=eval(ode,f(x) = x-1/2);
dsolve({ode1,y(0)=0,y(1)=0});

You will notice that now there are infinitely many solutions.

You have a polynomial system with 7 polynomials in the variables indets(sys,name) = {omega, x, a[-1], a[0], a[1], b[-1], b[0]}.

The coefficients are just integers. So in principle you can just try
solve(sys);
that will be understood by Maple as solve(sys,{omega, x, a[-1], a[0], a[1], b[-1], b[0]} );
Does the computation ever come to an end?
Try fsolve the same way, but be aware that you are only getting one solution out of many
fsolve(sys);
fsolve(sys, complex);

## You could also try SolveTools:-PolynomialSystem:

SolveTools:-PolynomialSystem(sys,{omega, x, a[-1], a[0], a[1], b[-1], b[0]},maxsols=20);

The following may not do what you intended, but at least the syntactical problems are gone.

You don't need any assumption on a.
 

restart;
N:=10;
b[0]:=1: b[1]:=-1: b[2]:=3: #Example
####
for n from 3 to N do    
  if (n mod 3) = 0  then
    b[n] :=b[0]   
  elif   (n mod 3) = 1 then
    b[n]:=b[1]   
  elif  (n mod 3) = 2  then
    b[n]:=b[2]   
  end if  
end do;
####
seq(b[i],i=1..10);
####
u[1]:=1: #Example
####
## You may want to assign a value to a before this, but you don't have to:
for n from 1 to N do  u[n+1]:=expand(a*u[n])+b[n]  end do;
### Seeing what you have for a = 0.4:
seq(eval(u,a=0.4)[n],n=1..N);

Beginning with Maple 2018, you can just do diff(A, x).
In earlier versions, but later than Maple 12 you can use diff elementwise: diff~(A,x).
In Maple 12 and earlier you can use map: map(diff, A, x).

map works in all versions. Elementwise works in all versions from Maple 13 and up.
Example:
 

A:=Matrix( [[sin(x),4*exp(3*x)],[cos(5*x),sinh(cos(x))]] );
diff(A,x);        # Maple 2018 and later
diff~(A,x);      # Maple 13 and later
map(diff,A,x); # All versions

The error you get in (e.g.) Maple 2017 from trying diff(A,x) is

Error, non-algebraic expressions cannot be differentiated

The error you get in Maple 12 from doing diff~(A,x) is more obscure:

Error, missing operator or `;`

Use value. It turns the inert Int into int.

Example:
 

V:=Vector[row]([seq(Int(x^n/L^(n+2),x=0..L),n=1..5)]);
value(V);

You are using evalm. That was used in the old and deprecated linalg package.

To get a matrix similar to yours you can do
 

V^+.V;
value(%);

Since you are using Maple 13, I tried in Maple 12 (don't have 13):
My guess is that there is not much difference between the two.
 

## Maple 12.
restart;
eqs:={(13/4)*m-(7/4)*n-3 = 0,-(17/2)*n*2^n +34*m= 0};
plots:-implicitplot(eqs,m=-1..3,n=-4..3);
sol:=solve(eqs);
evalf(sol); # Your desired solution, but as floats.
A:=allvalues(sol);
evalf(A); # Misses the other!
E:=eliminate(eqs,m);
eq_n:=op(E[2]);
plot(eq_n,n=-3..3); # Looks like 2 solutions (at least)
diff(eq_n,n);
z_n:=solve(%,n); 
#The derivative has one zero only. 
#Thus there cannot be more than 2 (real) solutions for n (and hence for (m,n) ).
evalf(z_n);
plot(diff(eq_n,n),n=-1..1);
fsolve(eq_n);
fsolve(eq_n,n=1);

So there are two real solutions for (m,n).