Your main problem is that you need to change linalg[matrix] to Matrix. Also, Matrix multiplication is with `.`, not `*`. Then you need to convert the Matrices to sets before passing them to fsolve. I also made some massive simplifications to your code that was not correcting errors.

restart:
Digits:= 20:

N:= 3:
 
y11:= x-> sum(a[n]*x^n, n= 0..N):
y12:= x-> sum(b[n]*x^n, n= 0..N):
y21:= x-> sum(c[n]*x^n, n= 0..N):
y22:= x-> sum(d[n]*x^n, n= 0..N):

A:= Matrix(2, 2, [1, -1, 1, exp(x)]):
B:= Matrix(2, 2, [-3*exp(-x)-1, 2-2*exp(-x), -3*exp(-x)-2, 1-2*cosh(x)]):
Z:= Matrix(2, 2, [3,0,1,1]):

Y:= Matrix(2, 2, [y11, y12, y21, y22](x)):
YY:= diff~(Y,x):

S:= `union`(seq(convert(eval((YY-A).(Y-B), x= n/N), set), n= 1..N)):
S1:= eval(S, convert(eval(Y, x= 0)=~ Z, set)):
sol:= fsolve(S1);

alias(h=f(x,y));

I think that Joe's abbreviated code is brilliant, but it still has time complexity O(|L|*|E|). Here's a version with time complexity O(|L|+|E|) which relies on Maple's table lookups using essentially constant time.

list_minus:= proc(L::list, E::{list,set})
local R:= table(sparse), e;
     for e in E do R[e]:= R[e]+1 end do;
     remove(e-> if R[e] > 0 then R[e]:= R[e]-1; true else false end if, L)
end proc:

Test:

list_minus([1,2,4,6,2,1,3,6,2], [7,4,2,5,2]);

Make the plot's color depend on the animation parameter. For example,

plots:-animate(plot, [x+b, x= -1..1, color= COLOUR(HUE, (b+1)/2)], b= -1..1, trace= 5);

See ?Optimization and ?Optimization,NLPSolve .

This is almost identical to your last question. The answer is identical: You need to use value to force evaluation of the integrals (and hence the removal of their dummy variable _z1). Change the line

f(tau):=sum((f[j])(tau),j=0..n);

to

f(tau):= value(sum(f[j](tau), j= 0..n));

My favorite way to plot a complex expression over a real parameter range is

plot([Re,Im](-x+sqrt(-5*x)), x= -10..10);

There's also ?plots,complexplot and ?plots,complexplot3d .

That form is called prefix form. In Maple 18, this is accomplished by

InertForm:-Parse("x*(y+z)/x");

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.