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

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

@acer Thanks for the explanation.
A comment:
map(eval,M+K);
does seem to evaluate all elements.

@acer Thanks for the explanation.
A comment:
map(eval,M+K);
does seem to evaluate all elements.

@Divhhh The time is given by the line integral

where s is arclength and v is velocity. From conservation of energy we have
1/2*m*v^2+m*g*y = C = 0 (if v = 0 at y = 0).
Thus v = sqrt(-2*g*y).
If the curve is given parametrically by x = x(u), y = y(u) (avoiding the use of t in order not to get confused), then ds/du = sqrt(x'(u)^2+y'(u)^2). Thus time spent on the curve (with u = 0 .. u0) is given by:

T = int( sqrt(x'(u)^2+y'(u)^2)/sqrt(-2*g*y(u)), u = 0..u0).
In Maple:

restart;
a := 1; b := -1:
g := 9.81:
eq1 := r*(u-sin(u)) = a; eq2 := r*(cos(u)-1) = b;
sol := fsolve({eq1, eq2}, {r, u});
r0,u0:=op(subs(sol,[r,u]));
x := u->r0*(u-sin(u));
y := u->r0*(cos(u)-1);
plot([x(u),y(u),u=0..u0]);
Int(sqrt(diff(x(u),u)^2+diff(y(u),u)^2)/sqrt(-2*g*y(u)),u=0..u0);
evalf(%);

@Divhhh The time is given by the line integral

where s is arclength and v is velocity. From conservation of energy we have
1/2*m*v^2+m*g*y = C = 0 (if v = 0 at y = 0).
Thus v = sqrt(-2*g*y).
If the curve is given parametrically by x = x(u), y = y(u) (avoiding the use of t in order not to get confused), then ds/du = sqrt(x'(u)^2+y'(u)^2). Thus time spent on the curve (with u = 0 .. u0) is given by:

T = int( sqrt(x'(u)^2+y'(u)^2)/sqrt(-2*g*y(u)), u = 0..u0).
In Maple:

restart;
a := 1; b := -1:
g := 9.81:
eq1 := r*(u-sin(u)) = a; eq2 := r*(cos(u)-1) = b;
sol := fsolve({eq1, eq2}, {r, u});
r0,u0:=op(subs(sol,[r,u]));
x := u->r0*(u-sin(u));
y := u->r0*(cos(u)-1);
plot([x(u),y(u),u=0..u0]);
Int(sqrt(diff(x(u),u)^2+diff(y(u),u)^2)/sqrt(-2*g*y(u)),u=0..u0);
evalf(%);

A change of variable works for both:

G:=Int(x*exp(-x)*exp(-exp(-x)), x = -infinity .. infinity);
IntegrationTools:-Change(G,{x=-ln(u)});
value(%);
V:=Int((x-gamma)^2*exp(-x)*exp(-exp(-x)), x = -infinity .. infinity);
IntegrationTools:-Change(V,{x=-ln(u)});
value(%);