For the first plot which doesn't work, the expression needs to be simplified to clear the units. (You seem to have already noticed that the output of breakforce needs simplification.) This can be corrected by changing the definition of breakforce from

breakforce:= unapply(...);

to

breakforce:= simplify@unapply(...);

Then you need to change the two xs in the plotted expression to x*Unit(MPa). That is, you need to change TS1(x) to TS1(x*Unit(MPa)) and change por1(x) to por1(x*Unit(MPa)).

For the second plot, you also need to change the xs as described above. And likewise for the next plot, the one of por1(xx).

I will answer your other questions, such as the one about displaying arrays of plots, in a separate Answer, when I get to it.

Extracting a coefficient from an integrand and placing it in front of the integral is done by the IntegrationTools:-Expand command.

J:= Int(a*f(x),x);

IntegrationTools:-Expand(J);

Are you sure that you're getting 10 digits? I think that you're getting 15 digits. Count them, please. By default, Matrix computations are done in hardware double precision if Digits <= 15. You change that default by setting UseHardwareFloats.

Set

UseHardwareFloats:= false:  Digits:= 5:

Why are you using MatrixPower? If you want A^11, just type A^11.

Assuming that n is a specific known integer, then

P:= (j,n)-> m*(m/2)^(j-1)/(2*(m/2)^n-1):
seq(P(j,n), j= 1..n-1);

By default Heaviside(0) is undefined. That's where your undefineds are coming from. You can change it to 1 by setting the environment variable _EnvUseHeavisideAsUnitStep. Also, the assumptions, combine, and simplify are not doing any good in this case.

restart;
_EnvUseHeavisideAsUnitStep:= true:
w123 := c1*r*Heaviside(r-a)+c2*r:
int(r*(int((diff(w123, r))^2/r, r)), r):

collect(%, [Heaviside, ln, r]):
convert(%, piecewise, r);

Note that subticks is not listed as one of the options in the error message. The subticks option can only be specified when the tickmarks option is part of an axis option. Like this:

plot(
     sin(x), x= -1..1, color= "DarkBlue",
     axis[1,2]= [tickmarks= [spacing(0.25), subticks= 2]]
);

Note also that I changed the color option. This color must be in quotes, and must be CamelCased.

If you include the option axes= frame in both plots and then rotate them so that you can see the scale of the z-axis, then I think that you'll understand what's happening. Every frame of the animation necessarily has the same z-axis, approximately -35..1. The static plot has z-axis approximately -0.6..0.3.  Let me know if you need more explanation.

Mehdi's animation worksheet had some formatting problems. Here is a different technique to get an animated W(x,t). In this worksheet, I use all your original variable names (indeed, I simply cut-and-pasted from your two original worksheets).

Download animated_beam.mw

To get an analytic colution with numeric coefficients, use the option convert_to_exact= false. Using Mehdi's formulation of the system, the command is

ans:= dsolve(sys union BCs, {X(t), Y(t)}, convert_to_exact= false);

This will give a solution that is 12 screens long, so there will be significant round-off error in the numeric evaluation. Indeed, if you want to evaluate it numerically, it would probably be better to use Mehdi's dsolve(..., numeric) solution. If you do evaluate the analytic solution at specific real values of t, you will get small imaginary parts. These are due to round-off error, and they should be ignored or discarded. Like this:

evalf(eval(ans, t= 1));

simplify(fnormal(%), zero);

You have five boundary conditions, but you should have only four.

It can be done like this, for example. This code does the same thing as the code that you posted above.

restart:
N := 10; #number of oscillators
K := 100; #number of times process is run
phi := 5; #bump given to oscillators when one fires
alpha := .5; #constant that affects extent to which size of firing populations affects bump
R:= rand(0..100); #100 should be a named constant
S := Array([seq(Record[packed]('val'= R(), 'sensitivity'= `?`, 'strength'= `?`), j= 1..N)]);

to K do
     t:= 100 - max(seq(S[j][val], j= 1..N));

     mu:= 0;
     for j from 1 to N do
          S[j][val]:= S[j][val] + t + phi;
          if S[j][val] > 99 then S[j][val]:= 0 end if;
          if S[j][val] = 0 then mu:= mu+1 end if
     end do;

     for j from 1 to N do
          if S[j]<>0 then
               S[j][val]:= S[j][val] + exp(alpha*mu);
               if S[j][val] > 99 then S[j][val]:= 0 end if
          end if
     end do
end do;

When you want to add a finite specific number of terms (rather than perform a symbolic summation) you should use add instead of sum. Then you will get the correct answer. By using sum, you are causing the premature evaluation of 0^(2*k), which evaluates to zero,  k being initially unknown. With add, the terms are not evaluated until the values of k are substituted, and then 0^0 evaluates to 1.

This issue has been discussed on Maple forums for decades.

Like your factor issue of today, there is also a workaround for this problem using frontend. This time we do freeze I, which frontend does by default. To freeze it means to treat it as just a name. In the factor issue, we froze the Pi (frontend froze it by default) but not the I.

frontend(convert, [[x, [I*x-1, 1]], parfrac, x], [{`+`,`*`,list}]);

Here is how to work around this issue with frontend:

frontend(factor, [Pi*(x^2+1), {I}], [{`+`, `*`}, {I}]);

I can't come up with logical argument about why things are the way they are with factor, except to say that perhaps it was too difficult for it to be designed otherwise, and that the designers knew that these idiosyncracies could be circumvented with frontend.

f5:= (a,b,c,d,e)-> a^2+b^2+c^2+d^2+e^2:
f3:= (ab, cd, e)-> ab[1]^2 + ab[2]^2 + cd[1]^2 + cd[2]^2 + e^2:
f5(1,2,3,4,5);
                                                   55
f3([1,2], [3,4], 5);
                                      55

There is a Maple procedure to automate this: codegen:-packparams:

codegen:-packparams(f5, [a,b,c], ABC);