Do not use a variable and an indexed variable with the same name, so I replaced  xi  with  xi_ :

restart;
l := -2;                             
m := 1;                             
k := sqrt(-1/(6*beta))/l;                      
w := (1/5)*alpha/(beta*l);                          
a[2] := -12*sqrt(-1/(6*beta))*alpha*m^2/(5*l*l);                   
a[0] := 0;                              
a[1] := 0;                              
F := -l*C[1]/(m*(C[1]+cosh(l*(xi_+xi[0]))-sinh(l*(xi_+xi[0]))));          
beta := -2;                              
alpha := 3;                               
C[1] := -1/1000;
xi[0] := 1;                               
xi_ := k*x-t*w;
u := a[0]+a[1]*F+a[2]*F*F;  
plot3d(u, x = -3 .. 3, t = -3 .. 3);

Do not use a variable and an indexed variable with the same name, so I replaced  f  with  F . Also made a few more fixes. Now there are no errors:

Download ADMTRY1_new.mw

Just use the  convert(x, float)  command or  evalf  command  for the finish result.

Use  CodeTools:-Usage  command for this. For ease of comparison, I increased the size of the array:
 

Download procedure_new1.mw

I think the problem can be solved in different ways. Below is an example of visualization of coin tosses. CoinToss procedure returns the result of N coin tosses in the form of a list of coordinates of the vertices of the broken line. Each segment up means the fall of one side of a coin, down - the other one. As an example, a simulation of 20 tosses of 2 independent coins is given. Animation of this is also shown.

Download CoinToss.mw

I think that it will be interesting for you to get the law of distribution of the number of points when throwing  n dices. In other words, it is necessary to find the probability that when throwing  n dice the sum of points will be equal to  S  for any possible S . For this we can use the procedure  Composition . Details on the use of this procedure can be found here.

Composition := proc (n::nonnegint, k::posint, res::{range, nonnegint} := 0)
local a, b, It, L0; 
if res::nonnegint then a := res; b := n-(k-1)*a  else a := lhs(res); b := rhs(res) fi;
if b < a or b*k < n then return `No solutions` fi; 
It := proc (L)
local m, j, P, R, i, N;
m := nops(L[1]); j := k-m; N := 0;
for i to nops(L) do
R := n-`+`(op(L[i]));
if R <= b*j and a*j <= R then N := N+1;
P[N] := [seq([op(L[i]), s], s = max(a, R-b*(j-1)) .. min(R, b))] fi;
od;
[seq(op(P[s]), s = 1 .. N)];
end proc;
L0 := [[]];
(It@@k)(L0); 
end proc:

Example of use - find all the ways to get a total of 12 points when throwing 3 dice:

Composition(12, 3, 1..6);
nops(%);      

   [[1, 5, 6], [1, 6, 5], [2, 4, 6], [2, 5, 5], [2, 6, 4], [3, 3, 6], [3, 4, 5], [3, 5, 4], [3, 6, 3], [4, 2, 6], [4, 3, 5], [4, 4, 4], [4, 5, 3], [4, 6, 2], [5, 1, 6], [5, 2, 5], [5, 3, 4], [5, 4, 3], [5, 5, 2], [5, 6, 1], [6, 1, 5], [6, 2, 4], [6, 3, 3], [6, 4, 2], [6, 5, 1]]
                                                                     25

In the example below we obtain the law of distribution of the amount of points when rolling 4 dice :

n:=4:
# Probabilities
P:=[seq([S,nops(Composition(S,n,1..6))/6^n], S=n..6*n)];
# Check
add(P[k,2], k=1..nops(P));

Edit.

Your integral equation has the exact solution  u(x)=1+3/4*x . The solution using the trapezoid rule for  N=100  gives a pretty good approximation:

Download IntEq.mw

I did not find a way to increase the size of this font inside Maple. As a workaround, you can use a Windows screen magnifier or just wear glasses. See the screenshot:

If you have specific points, you can use  geom3d  package.

Example:

restart;
with(geom3d):
point(P, 5, 6, 7);
point(A, 1, -2, 3);
point(B, 2, 4, 1);
line(AB, [A,B]);
projection(Q, P, AB);
detail(Q);
                               

A:=isolve(n+k=2*c);
is(A<>NULL and eval(n, A)::name);
                          A := {c = _Z1, k = 2*_Z1-_Z2, n = _Z2}
                                                       true

Above  _Z1  and  _Z2  are any integers.

You can use  PolyhedralSets  package to get the vertices and edges of your polygon. This also works in the multidimensional case.

Your example (failed to load the contents of the file, so the image is shown, the link to the file below):

You can extract all vertices in one list with  Relations  command:
V:=Vertices(ps):
Relations~(V);
map(v->eval([x,y],v), %);
                  [[y = -1, x = -1], [y = 1, x = -1], [y = -1, x = 1], [y = 1, x = 1]]
                                       [[-1, -1], [-1, 1], [1, -1], [1, 1]]

The same for edges.

Download PS.mw

Use  factors  command for this:

f:=x*(x+1)^4*(x^2+x+1)*(x^3+x^2+1)^5:
factors(f);
map(t->op(1,t), %[2]);
sort(%);
                         

I usually use  is  command for such checks:

is(2*cos(phi)^2-1=cos(2*phi));
                                                             true

And more often I look  here , if I don’t know any rare identity.

Edit.

Maple 2018.2 returns the result in the form of finite sums and products (without hypergeometric functions):

restart;
int(sin(x)^n, x) assuming n::posint;
eval(%, n=10);
value(%);
int(sin(x)^10, x); # The direct calculation

restart;
Ke := Matrix(4, 4, {(1, 1) = (12*I)*E/l^3, (1, 2) = (6*I)*E/l^2, (1, 3) = -(12*I)*E/l^3, (1, 4) = (6*I)*E/l^2, (2, 1) = (6*I)*E/l^2, (2, 2) = (4*I)*E/l, (2, 3) = -(6*I)*E/l^2, (2, 4) = (2*I)*E/l, (3, 1) = -(12*I)*E/l^3, (3, 2) = -(6*I)*E/l^2, (3, 3) = (12*I)*E/l^3, (3, 4) = -(6*I)*E/l^2, (4, 1) = (6*I)*E/l^2, (4, 2) = (2*I)*E/l, (4, 3) = -(6*I)*E/l^2, (4, 4) = (4*I)*E/l}):
k:=E*I/l^3:
``(k).(Ke/k);