I suppose the problem Maple has with this example is the error handling?

The integrand G

G:=diff(-ln(1-Fy), y)*fy;
# is clearly a smooth function defined on all of R and falls of very rapidly as |y| -> infinity
plot(G,y=-10..10);
asympt(G,y);
#I got the following slightly disturbing error from this command (Digits = 10):
evalf(Int(G,y=-10..10));
# error message:
Error, (in property/LinearProp/+) too many levels of recursion
#However the results of the following were OK:
evalf[15](Int(G, y = -5 .. 5));
evalf[15](Int(G, y = -10 .. 10));
evalf[15](Int(G, y = -20 .. 20));

That weird output makes me ask: What was the input?
Here is an example:

with(LinearAlgebra):
A:=Matrix([[5, -3, 3], [-4, 8, -7], [-5, 3, 2]]); #Simple example
Eigenvalues(A);
evalf(%);

Try your loop on
T:={A,B,C};
and you will see that you are missing the union of all 3 sets A, B, and C.

Try your loop on
T:={A,B,C};
and you will see that you are missing the union of all 3 sets A, B, and C.

The answer from the following is not quite satisfactory:

completesq(ex,[alpha-delta,beta-delta]);

since the sum of the last 3 terms is zero. This could be handled I'm sure, but the main problem remains: How to find a good second argument for completesq (without knowing the final answer in advance)?
I suspect that there is no general answer to that question.

Note:

The unsatisfactory behavior mentioned really is due to Student:-Precalculus:-CompleteSquare:

u:=expand(  (x-a+b)^2+(y-c)^2);
Student:-Precalculus:-CompleteSquare(u,[x,y]);

but could in the case where Student:-Precalculus:-CompleteSquare is called with more than one argument be handled by expanding (or simplifying) the terms not involving the arguments 2..-1 as illustrated in this version of your procedure:

completesq:=proc(ex::algebraic,l::{algebraic,list(algebraic)})
 local l1:=freeze~(l),
 s:=zip(`=`,l,l1),
 s1:=foldr(algsubs,expand(ex),`if`(s::list,op(s),s)),
 s2:=Student:-Precalculus:-CompleteSquare(s1,l1);
 s2:=evalindets(s2,freeof(l1),expand);#Using expand
 thaw(s2);
end proc:

The answer from the following is not quite satisfactory:

completesq(ex,[alpha-delta,beta-delta]);

since the sum of the last 3 terms is zero. This could be handled I'm sure, but the main problem remains: How to find a good second argument for completesq (without knowing the final answer in advance)?
I suspect that there is no general answer to that question.

Note:

The unsatisfactory behavior mentioned really is due to Student:-Precalculus:-CompleteSquare:

u:=expand(  (x-a+b)^2+(y-c)^2);
Student:-Precalculus:-CompleteSquare(u,[x,y]);

but could in the case where Student:-Precalculus:-CompleteSquare is called with more than one argument be handled by expanding (or simplifying) the terms not involving the arguments 2..-1 as illustrated in this version of your procedure:

completesq:=proc(ex::algebraic,l::{algebraic,list(algebraic)})
 local l1:=freeze~(l),
 s:=zip(`=`,l,l1),
 s1:=foldr(algsubs,expand(ex),`if`(s::list,op(s),s)),
 s2:=Student:-Precalculus:-CompleteSquare(s1,l1);
 s2:=evalindets(s2,freeof(l1),expand);#Using expand
 thaw(s2);
end proc:

@mehdi_mech If the realistic model really behaves like the one you present, then a lot of oscillation takes place over a very short time (t).
The following much simpler example shows the problem:

restart;
eq:=diff(a(t),t,t)+10^12*a(t)=0;
res:=dsolve({eq,a(0)=0,D(a)(0)=1}); #Analytical solution available here
plot(rhs(res),t=0..1e-5);
#You see the problem for numerical computations if you change the t-interval above to e.g. t=0..1
#Numerical solution:
resnum:=dsolve({eq,a(0)=0,D(a)(0)=1},numeric);
resnum(1e-5); #OK with a tiny t-value, but try a bigger one.
#As in my answer above we could make a change of time variable, but that won't make it easier to get values for a(t) for much larger values of t. The point in scaling time was mainly to show that quite a lot goes on over a very short time, so solving over a comparatively long interval is bound to take a long time.
In fact once you realize this, you may as well not scale time.

@mehdi_mech If the realistic model really behaves like the one you present, then a lot of oscillation takes place over a very short time (t).
The following much simpler example shows the problem:

restart;
eq:=diff(a(t),t,t)+10^12*a(t)=0;
res:=dsolve({eq,a(0)=0,D(a)(0)=1}); #Analytical solution available here
plot(rhs(res),t=0..1e-5);
#You see the problem for numerical computations if you change the t-interval above to e.g. t=0..1
#Numerical solution:
resnum:=dsolve({eq,a(0)=0,D(a)(0)=1},numeric);
resnum(1e-5); #OK with a tiny t-value, but try a bigger one.
#As in my answer above we could make a change of time variable, but that won't make it easier to get values for a(t) for much larger values of t. The point in scaling time was mainly to show that quite a lot goes on over a very short time, so solving over a comparatively long interval is bound to take a long time.
In fact once you realize this, you may as well not scale time.

Matrices do not in general commute. So even if both products AB and BA exist then they are unlikely to be equal.
Now if A.X = B and if X has an inverse then A = B.X^(-1) (and not X^(-1).B ) .
A reason for not using the division sign / in the context of matrices is: Which of  B.X^(-1) or X^(-1).B should B/X be? This is much more of a problem if / is replaced by the horizontal version:

Should this mean A.X^(-1) or X^(-1).A ?

Matrices do not in general commute. So even if both products AB and BA exist then they are unlikely to be equal.
Now if A.X = B and if X has an inverse then A = B.X^(-1) (and not X^(-1).B ) .
A reason for not using the division sign / in the context of matrices is: Which of  B.X^(-1) or X^(-1).B should B/X be? This is much more of a problem if / is replaced by the horizontal version:

Should this mean A.X^(-1) or X^(-1).A ?

@FKil w1 is a product of the positive constant F and a function of x, i.e. w1 = F*f(x). As I understand it you want to solve
max( F*f(x), x=0..l) = 1.43e6 w.r.t. F. Since F is constant and positive
max( F*f(x), x=0..l) = F*max( f(x), x=0..l).
Now max( f(x), x=0..l) can be found by Optimization:-Maximize as mentioned earlier. If the maximum is called r, then you have F*r = 1.43e6, so F = 1.43e6/r. Since r is found to be r = -0.237256300503585 we get F = -6.027237199*10^6.

@FKil w1 is a product of the positive constant F and a function of x, i.e. w1 = F*f(x). As I understand it you want to solve
max( F*f(x), x=0..l) = 1.43e6 w.r.t. F. Since F is constant and positive
max( F*f(x), x=0..l) = F*max( f(x), x=0..l).
Now max( f(x), x=0..l) can be found by Optimization:-Maximize as mentioned earlier. If the maximum is called r, then you have F*r = 1.43e6, so F = 1.43e6/r. Since r is found to be r = -0.237256300503585 we get F = -6.027237199*10^6.

@FKil In your worksheet I see the expression (which I shall call w):

w:=evalf(-N/(b*h)+((-EI*(diff(u(x), x, x)))*.5)*h/(0.10e-10*Iice));
#I have used evalf on it to make it shorter
#Collecting w.r.t. F
w1:=collect(w,F);
whattype(w1); #The result is a product
nops(w1); #Two factors
op(2,w1); #One factor is F and the other doesn't contain F:
indets(op(1,w1),name);

So maximizing op(1,w1) (same as maximizing w1 with F=1) gives you a real number (say r), and you want to solve r*F = 1.43e6.


So what is the problem?

@FKil In your worksheet I see the expression (which I shall call w):

w:=evalf(-N/(b*h)+((-EI*(diff(u(x), x, x)))*.5)*h/(0.10e-10*Iice));
#I have used evalf on it to make it shorter
#Collecting w.r.t. F
w1:=collect(w,F);
whattype(w1); #The result is a product
nops(w1); #Two factors
op(2,w1); #One factor is F and the other doesn't contain F:
indets(op(1,w1),name);

So maximizing op(1,w1) (same as maximizing w1 with F=1) gives you a real number (say r), and you want to solve r*F = 1.43e6.


So what is the problem?

@FKil Isn't it as simple as the following or am I misunderstanding something?

sm:=Optimization:-Maximize(eval(-N/(b*h)+((-EI*(diff(u(x), x, x)))*.5)*h/(0.10e-10*Iice), F = 1), x = 0 .. l);
eq:=sm[1]*F=sigma[t];
solve(eq,F);