When calculating the derivative of a function at a point, it is better to use the differentiation operator  D  than  eval  command (simpler and more compact syntax). For example:

eval(diff(x^2,x), x=2.5);
D(x->x^2)(2.5);
                                                 5.0
                                                 5.0

In your example, you are trying to differentiate the function at points where it is not defined. See:

restart;
y := x-> (1-Heaviside(x-0.10e-1))*(cos(9.0218219*x)-.99533285*sin(9.0218219*x)-.99999991*cosh(9.0218219*x)+.99533285*sinh(9.0218219*x))+(Heaviside(x-0.10e-1)-Heaviside(x-0.418e-1))*(.24369100*cos(7.7520047*x)-.36109859*sin(7.7520047*x)-.23739778*cosh(7.7520047*x)+.19615343*sinh(7.7520047*x))+Heaviside(x-0.418e-1)*(.95680995*cos(9.0218219*x)-.80884870*sin(9.0218219*x)-1.0381918*cosh(9.0218219*x)+1.1704059*sinh(9.0218219*x)):
convert(y(x), piecewise);
y(0.418e-1);
y(0.1e-1);

This bug occurs when using the palette in 2d math input. Instead, enter a vector or a matrix using angle brackets directly from the keyboard. It's much faster:

We can immediately (without even using  combinat  package) obtain all such sums for any k  (the list  L ) if we first construct the corresponding generating function (the polynom  P ).

Example:

n:=5:
P:=collect(expand(mul(x+~[a||(0..n-1)])), x);
L:=[seq(coeff(P,x,n-k), k=1..n)];
   

I guess that you are looking for real roots  (x1, y1) , and  t, b, n  are parameters. Your system has the form  {eq1, eq2} = {f(x1)=f(y1), g(x1)=g(y1)}  in which  f = x->t/(x-1)/x^2-b*n/6/x^2*x^(-n-1), g = x->t*(0.74/(x-1)/x^3-3*ln(x-1))+b/x^(n)*(1/2-0.74*n/6/x^4) . Obviously, this system always has an infinite number of solutions  (x1,y1) = (x,x) , in which x>1 . If  f  or  g  is monotonous, then there are no other solutions.

The  extrema  command allows in many cases to symbolically find local extremes of smooth multidimensional problems with smooth constraints. For the solution, the Lagrange multipliers method is used. We get a more accurate solution with less cost. Unfortunately, the exact symbolic result is quite cumbersome and its approximation is shown below:

Digits:=25:
min(extrema((x[1]-1)^2+(x[2]-1)^2+(x[3]-1)^2, x[1]*(x[2]^2+x[3]^2+1)=8, {x[1], x[2], x[3]}, 's')):
evalf(%);
evalf(s):
select(p->`and`(seq(eval(x[i], p)>=1 and eval(x[i], p)<=2, i=1..3)), %)[];
evalf[23](eval((x[1]-1)^2+(x[2]-1)^2+(x[3]-1)^2, %)); # Check

The output:
                                      0.69329028371527531623476
         {x[1] = 1.531665799394114675166437, x[2] = 1.453112437188548339089389, x[3] = 1.453112437188548339089389}
                                      0 .69329028371527531623477

From a geometrical point of view, we found the value of the objective function  (x[1]-1)^2+(x[2]-1)^2+(x[3]-1)^2  when its level surfaces  (x[1]-1)^2+(x[2]-1)^2+(x[3]-1)^2=С  (in your example these are ordinary spheres) touch the surface  x[1]*(x[2]^2+x[3]^2+1)=8  for the first time.

Addition. Since Optimization:-Minimize and extrema commands generally return local extremes (and extrema command - the critical points), it is useful to verify the result obtained using  DirectSearch:-GlobalOptima  command, which searches for the global extremum using a direct search without derivatives:

restart;
DirectSearch:-GlobalOptima((x[1]-1)^2+(x[2]-1)^2+(x[3]-1)^2, {x[1]*(x[2]^2+x[3]^2+1)>=8, seq(op([x[i]>=1,x[i]<=2]), i=1..3)});
      [0.693290285423826, [x[1] = 1.53164167847983, x[2] = 1.45310861919799, x[3] = 1.45314455787912], 411]       

Download CoinToss_new.mw

I don’t know how to extract the vector structures you are writing about, but just to simplify the trigonometric expression  applyrule command can help. First we apply  simplify command to your expression, and then applyrule command, and in the list we indicate the formulas that applyrule may need. The following is a simplified example of yours, in which I just replaced the names of the angles using simple indexing:

Expr:=cos(theta__a)*cos(theta__b)+sin(theta__a)*cos(phi__a)*sin(theta__b)*cos(phi__b)+sin(theta__a)*sin(phi__a)*sin(theta__b)*sin(phi__b):
# Simplification of Expr
applyrule([cos(t::anything)*cos(s::anything)+sin(t::anything)*sin(s::anything)=cos(t-s), cos(t::anything)*cos(s::anything)-sin(t::anything)*sin(s::anything)=cos(t+s), sin(t::anything)*cos(s::anything)+cos(t::anything)*sin(s::anything)=sin(t+s), sin(t::anything)*cos(s::anything)-cos(t::anything)*sin(s::anything)=sin(t-s)],simplify(Expr));

Output:       

In the code for  applyrule  I wrote down 4 well-known addition formulas for sine and cosine, but you can add formulas of multiple angles, etc.
The commands  applyop , subsindets  can also be useful for selective use.

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: