You have found a bug in Maple, and I'll give you a workaround in a moment. But your code itself also has 1 serious bug and 1 serious shortcoming, and it amazes me that these didn't confound the issue of the true bug that you found.

Bug: Change
f:= x-> piecewise(0 <= t, 1, t < 0, t - 1)
to
f:= t-> piecewise(0 <= t, 1, t < 0, t - 1)

Shortcoming: Change
F:= x-> int(f(t), t= 1..x)
to
F:= unapply(int(f(t), t= 1..x), [x])

Workaround: Change the unapply to
F:= unapply(int(f(t), t= 1..x, method= FTOC), [x])

So, why use unapply? Because F contains a symbolic operation, the integration, that can be performed independently (as well as once-and-for-all) of F's parameter x.

 

The Taylor or Maclaurin series for a function of a square matrix argument is the same as the series for the same function with a complex argument. The general non-commutativity of matrix multiplication plays no role in this because all powers of a single square matrix do commute. Maple can find the general term of the series for many fairly simple functions and express those series in sigma notation, like this:

convert(exp(A), FormalPowerSeries); #The second argument can be abbreviated FPS.

The commands series or taylor (essentially the same thing) can give you any finite number of terms of the series for nearly any function that can be represented as a power series, but they won't give the general term:

series(exp(A), A=0, 9);

Here is a procedure for computing the nth partial sum of the series for exp(A). The simplicity of this formulation amazes me! It avoids computing high powers of the matrix (which can easily lead to numeric overflow) by combining the updating of the matrix power and the updating of the factorial denominator with the single operation Ak.= A/k. This is an updating assignment operator, introduced in Maple 2018. It's equivalent to Ak:= Ak.(A/k).

#Computes the nth partial sum of the Maclaurin series for exp(A):
Exp:= (A::Matrix(square), n::nonnegint)->
local k, Ak:= rtable(identity, rtable_dims(A), rtable_options(A));
    Ak + add((Ak.= A/k), k= 1..n)
:    
#Examples:
Digits:= 32:
LA:= LinearAlgebra:
Display:= x-> (x, print(evalf[3]~(x)))[1]:
randomize():
A:= (Display@LA:-RandomMatrix)(
    (6,6), generator= rand(-2.0..2.0), datatype= sfloat
):
              [ 0.621    1.09     1.55  0.810  -0.464    0.198]
              [                                               ]
              [ -1.41   0.573  -0.0609  -1.66  0.0230    0.816]
              [                                               ]
              [ 0.689  -0.178     1.51  0.384    1.73  -0.0984]
              [                                               ]
              [  1.60  -0.448    0.815   1.79   -1.90    -1.53]
              [                                               ]
              [-0.209  -0.588    -1.19  -1.27   0.338     1.36]
              [                                               ]
              [ -1.73   0.935    0.441  -1.87   0.344  -0.0294]

E1:= Display(CodeTools:-Usage(Exp(A, 99))):
memory used=10.54MiB, alloc change=0 bytes, 
cpu time=62.00ms, real time=68.00ms, gc time=0ns

                [ 3.89     1.54   6.33   2.37   1.29  -0.195]
                [                                           ]
                [-20.6    0.996  -18.2  -23.9   7.80    11.1]
                [                                           ]
                [0.660   -0.262   3.31  -1.55   4.30    2.03]
                [                                           ]
                [ 30.8    0.288   29.1   37.0  -12.3   -17.5]
                [                                           ]
                [-14.5   -0.328  -14.5  -17.3   5.83    8.49]
                [                                           ]
                [-21.0  -0.0767  -18.4  -24.7   8.69    12.1]

#Compare:
E2:= Display(CodeTools:-Usage(LA:-MatrixExponential(A))):
memory used=2.85MiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=14.00ms, gc time=0ns

                [ 3.89     1.54   6.33   2.37   1.29  -0.195]
                [                                           ]
                [-20.6    0.996  -18.2  -23.9   7.80    11.1]
                [                                           ]
                [0.660   -0.262   3.31  -1.55   4.30    2.03]
                [                                           ]
                [ 30.8    0.288   29.1   37.0  -12.3   -17.5]
                [                                           ]
                [-14.5   -0.328  -14.5  -17.3   5.83    8.49]
                [                                           ]
                [-21.0  -0.0767  -18.4  -24.7   8.69    12.1]

#absolute and relative errors:
Display(<[Error__||(abs,rel)]=~ LA:-Norm~([E1-E2, (E1-E2)/~E2])>):
                          [                    -29]
                          [Error__abs = 1.72 10   ]
                          [                       ]
                          [                    -30]
                          [Error__rel = 3.99 10   ]

#My procedure Exp can handle symbolic matrices just as well. Here is MacDude's
#example done by it:
Exp(<a, b; c, d>, 2);
                [        1  2   1           1       1      ]
                [1 + a + - a  + - b c   b + - a b + - b d  ]
                [        2      2           2       2      ]
                [                                          ]
                [     1       1                1       1  2]
                [ c + - c a + - d c    1 + d + - b c + - d ]
                [     2       2                2       2   ]

If the square matrix A is diagonalizable (which is true for most numeric matrices in practice), then there's a better way than partial sums to compute f(A): If A = P^(-1).D.P where D is a diagonal matrix, then f(A) = P^(-1).f(D).P where f(D) can be computed simply by applying the complex-valued f to the entries on the diagonal. This is why the MatrixExponential command shown above is faster.

Many unexpected things can happen with alias because it acts at the level of what we call automatic simplification, which is one of Maple's first steps in the evaluation of an expression and often overlooks the intended meaning of an expression that would otherwise be obvious at the deeper levels of evaluation^[*1]. Thus, alias must be used with great care. It is often one of the first commands taught to a student just learning Maple, as it was for me. That is perhaps not a good idea. From years of answering Questions here, I've seen that problems often occur when trying to mix alias with indexed structures (either explicitly indexed such as a[1] or indexed-in-essence such as a1, a2)^[*1].

Here's an alternative to alias that still gives you all the benefits of abbreviation:

n:= 3:
A:= [a||(1..n)](r);
A[2];
dA:= diff(A,r);

Although it's not an issue in this case, one should also keep in mind that variable names constructed by concatenation are always global, regardless of whether the construction is done by ||, cat, nprintf, parse, or convert(..., name), regardless of whether the base name (a in this case) is local or global, and regardless of whether the construction occurs inside a procedure.

[*1] This paragraph from the help page ?alias explains a lot of this:

• Because aliases are resolved at the time of parsing the original input, they substitute literally, without regard to values of variables and expressions.  For example, alias(a[1]=exp(1)) followed by evalf(a[1]) will replace a[1] with exp(1) to give the result 2.718281828, but evalf(a[i]) will not substitute a[i] for its alias even when i=1.  This is because a[i] does not literally match a[1].
   

@Zeineb All 4 definite integrals from your most-recent worksheet compute very quickly to exact real values in Maple 2020. In your Maple, try doing it by limits, such as:

I3:= int(f*P3, x); 
simplify(eval(I3, x= 2) - limit(I3, x= 0, right));

This is also worth trying, and is less typing if it works:

int(f*P3, x= 0..2, method= FTOC); 

My guess is that the two methods above are equivalent. FTOC stands for Fundamental Theorem Of Calculus.

Most of the y-values have imaginary parts of significant magnitude.

I don't have Maple 2021 to check this, but I suspect that the problem that you report applies to any array plot. To test my hypothesis, try exporting this 2x2 array of plots to latex:

plots:-display(plot~(<x, x^2; x^3, x^4>));

I think that @mmcdara has made this seem more "complex" than it really is. Here is how to do the 3D plot that you want:

restart
:
par:= {
    phi= 0, lambda= 0.1, N= 5,
    M= sqrt(N*(N+1))*exp(I*phi), omegap= 10
}:
var:= {n,u,z}(t):
dsys:= {
    diff(n(t),t) =
        -2*(n(t)-N)+(u(t)-M)*exp(-2*I*omegap*t/lambda) + 
        ((z(t)-conjugate(M))*exp(2*I*omegap*t/lambda)),
    diff(u(t),t) = 
        -2*(1-I*delta)*u(t) + 2*(n(t)-N)*exp(2*I*omegap*t/lambda)+2*M,
    diff(z(t),t) = 
        -2*(1+I*delta)*z(t) +
        2*(n(t)-N)*exp(-2*I*omegap*t/lambda)+2*conjugate(M)
}:
res1:= dsolve(
    eval[recurse](dsys, par) union (eval(var, t= 0)=~0), 
    numeric, parameters= [delta], output= listprocedure
):

#Function to be plotted, and a helper function to speed things up:
Ndelta:= proc(delta::numeric)
option remember;
    res1(parameters= [delta]);
    #Return the instantiated procedure:
    sscanf(sprintf("%m", eval(eval(n(':-t'), res1))), "%m")[]
end proc
:
N:= (t, delta)-> 
    `if`([args]::list(numeric), Ndelta(delta)(t), 'procname'(args))
:
P3:= CodeTools:-Usage(
    plot3d(
        N(t, delta), t= 0..1, delta= -10..10, axes= boxed,
        tickmarks= [3, [-8, 8], 5],  
        titlefont= [Helvetica, roman, 18], 
        labels= [t, delta, n(t)], 
        labelfont= [Helvetica, roman, 24]
    )
);
memory used=8.47GiB, alloc change=0 bytes,
cpu time=84.97s, real time=78.92s, gc time=14.45s

#Show delta variation for fixed t:
P2:= CodeTools:-Usage(
    plot(N(.5, delta), delta= -10..10, numpoints= 49, adaptive= false)
);
memory used=4.69GiB, alloc change=48.00MiB, 
cpu time=50.72s, real time=47.10s, gc time=8.64s

At this scale, the variation with respect to delta is visually insignificant compared to the variation with respect to t. You can see it by using a much smaller t range or, perhaps, a much larger delta range.

I'm not sure what you mean by "system of equation". What you show is a single piecewise expression in two variables. You got an error because there are two independent free variables. In this case, they can be easily replaced by a single variable. So, perhaps the below plot is what you want. Let me know.

pw:= piecewise(t-tau <= 0, 2, 2+5*(t-tau)):
plot(simplify(pw, {t-tau = s}), s= -2..2, 0..10, labels= [t-tau, ``]);

I'm not pleased about the default font choice for the tau on the axis label! I suspect that I could easily fix that if you want. Also, it may appear as a normal Greek tau on your display (not here, but within Maple).

Here's an example of what @acer meant by "rescaling and forced tickmarks". For use over significantly wider or narrower ranges, the round in procedure ReTick will need to be changed. [Edit: I now see that acer has updated his Answer to provide such an example. His acts after the plot is created, and mine acts before. Both approaches have their pros and cons. In either case, having the z-axis constrained to a "nice" range by view makes things significantly easier.]

ReTick:= proc(ab::range(realcons), s::positive, n::posint:= 5)
local a:= lhs(ab), h:= (rhs(ab)-a)/n/2, k;
    ':-tickmarks'= 
        [seq](round(a+k*h)/s = `if`(k::odd, "", round(a+k*h)), k= 0..2*n)
end proc
:
(sx,sy):= (1,5): #axes scale factors
f:= (x,y)-> abs(Zeta(x+y*I)):
(xr,yr):= (-4..5, -10..40): #axes ranges
plot3d(
    f@((x,y)-> (sx*x, sy*y)), map(`/`, xr, sx), map(`/`, yr, sy),
    labels= ['x', 'y', ``], view= [DEFAULT, DEFAULT, 0..4],
    axis[1]= [ReTick(xr, sx, 4)], axis[2]= [ReTick(yr, sy)],
    orientation= [-15, 75, 0], scaling= constrained
);

In Maple's 2D Input (the default input mode, which you show in your Question):

1. When operands are immediately adjacent (with no space between them), there is no implied multiplication. If the second operand is in parentheses (as in your case), Maple interprets it as a function argument for which the first operand is the function. This interpretation is likely to be syntactically correct (because Maple is willing to use almost anything as a function (see help page ?evalapply)), so, unfortunately, it often doesn't produce an error message. This is likely the number-one cause of errors related to the input mode.
2. When there is a space or spaces, an implicit multiplication is automatically inserted. Sometimes, it becomes visible; whether or not it's visible isn't of much importance other than for superficial debugging. For any debugging requiring more than a few minutes, you'd convert the input to 1D anyway.
3. Explicit multiplication operators can also be used.

In 1D input (aka Maple Input or Maple Notation), as shown by Kitonum:

1. There is never implicit multiplication; all multiplication is through explicit operators.
2. When operands are adjacent with or without space, and the second operand is in parentheses, the function/argument interpretation described in the last section will be used.
3. The code can be read, understood, and debugged by any human reader who knows Maple. There are no subtle questions of interpretation.

The following short code will convert your L into a Maple matrix of your specified form:

L:= [[[1,2],3], [[1,3],4], [[1,4],4], [[1,7],7]]:
ExMat:= <<` A ` |` B `>, <<(``@op)~(L[.., 1])> | <L[.., 2]>>>;

I don't know about the ability of Excel to accept this form. It may be that conversion to strings will be needed, as Tom did. If so, that conversion can be just as easily done to the ExMat created above. I wrote this Answer primarily because the above ExMat may be sufficient for you to do whatever you want to do directly in Maple.

Code notes: 

1. ``@op is a composed operator (with @ being the composition operator (see ?@)) that converts a list into an equivalent form with "hard" parentheses. Technically, this form is called a function by Maple, with the function's name being the empty name ``. The commands expand or op will return the underlying sequence, which'll also make the parentheses disappear. Unlike in a list, in a vector or matrix a "naked" sequence (i.e., one without parentheses) can be a single entry. If this is fine for your purpose, you can simplify (``@op) to just op, which is the most-basic deconstructor of Maple objects (and likely the most commonly used Maple command with an alphabetic name).
2. <... | ... | ... > is a constructor of row vectors, and also a columnwise constructor for matrices when its arguments all have the same number of rows.
3. < ..., ..., ... > is a constructor of column vectors, and also a rowwise constructor for matrices when its argument all have the same number of columns.
4. If there's ambiguity between operators 2 and 3 because there's no |s and no commas, then it's 3. If LL is a list (possibly containing sublists) then <LL> forms the column vector of the entries. So, this implicitly removes the outer (only) level of [ ] that make LL a list.
5. Operators 2 and 3 can be arbitrarily nested; I believe that the result is always either a vector or a matrix, and that these operators cannot be used to build higher-dimensional arrays.
6. ~ is the elementwise operator (see help page ?elementwise). 
7. Just like a matrix, a list of lists can be multi-indexed as shown: L[.., 1]; however, if this were done on a large scale, it's not as efficient as it is for a matrix.
8. When the "naked" (i.e., with no apparent operands) range operator .. is used as an index, it refers to the entirety of the dimension corresponding to its position. In some other languages (Matlab, Python), a colon is used for this purpose. I think that this is the only Maple infix operator that can be used "naked". Technically, when it's used like this, the operands are NULL, which can also be represented ().
9. `...` is the symbol-forming operator, which is used to include non-alphanumeric characters in a symbol. In this case, I used it to pad the A and B with spaces. This serves two purposes (in this case): 1) display formatting, and 2) it protects against the possibility that A or B have been assigned values. Maple symbols are similar to strings, the main differences being 1) symbols can be assigned values, and 2) the two types may display differently, although the exact nature of that difference depends on where they're being displayed. 

Perhaps you should rethink your data structures. Why use an Array in the first place? I can't say for sure without knowing more details of your application, but from all that you've said so far, I think that a table would work better for you. 

If you change 0.5*X1 + 0.5*X2 to (X1+X2)/2, then it'll work without needing any change to Digits. This'll also work for Beta(5,5) and Beta(10,10).

Here is a procedure for it. The procedure is a tour de force of Maple structured-type programming, which is a somewhat esoteric topic, so it may be difficult to understand. The main idea is that by converting derivatives to D form, the order is easier to determine.

specdifforder:= proc(
    e, o::type, u::{name, function, set({name, function})}:= {}
)
local
    O:= {`if`(o::set, o, {o})[], identical(D)},
    U:= identical(
        `if`(
            u={},
            indets(e, function(name)),
            ((S,R)-> indets(e, specfunc(S)) union R)
                (selectremove(type, `if`(u::set, u, {u}), name))
        )[]
    ),
    Op:= ()-> curry(op, args),
    uT0:= And(
        U,    
        Not({specfunc(D), specop(`@@`) &under Op(0)}) &under Op(0)
    ),        
    uT:= And(
        U &under (f-> op([0, ..], f)(op(f))),
        Or([O], O &under nops) &under [Op([0, 0, ..])]
    ),
    R:= {}      
;
    if 0::O then R:= indets(convert(e, D), uT0) fi;
    R union indets(e, And(specfunc(diff), uT &under (convert, D)))
end proc
:
#Examples
e:= 
     diff(u(x,y),x,y) + diff(u(r,s),r,r) + u(x,y) + diff(u(x,y),x) + 
     diff(u(x),x) + diff(u(x),x,x)
:
specdifforder(e, 0, u);
                                  {u(x, y)}
specdifforder(e, 0);
                                  {u(x, y)}
specdifforder(e, 1, u);
                    / d         d            d         \ 
                   { --- u(x), --- u(r, s), --- u(x, y) }
                    \ dx        dr           dx        / 

specdifforder(e, 1, u(r,s));
                                / d         \ 
                               { --- u(r, s) }
                                \ dr        / 
specdifforder(e, {0,1}, u);
                / d         d            d                  \ 
               { --- u(x), --- u(r, s), --- u(x, y), u(x, y) }
                \ dx        dr           dx                 / 

specdifforder(e, 2);
                  /  2          2              2          \ 
                  | d          d              d           | 
                 < ---- u(x), ---- u(r, s), ------ u(x, y) >
                  |   2          2           dy dx        | 
                  \ dx         dr                         / 

specdifforder(e, 2, u(x));
                                 /  2      \ 
                                 | d       | 
                                < ---- u(x) >
                                 |   2     | 
                                 \ dx      / 

specdifforder(e, 2, {u(x), u(r,s)});
                          /  2          2         \ 
                          | d          d          | 
                         < ---- u(x), ---- u(r, s) >
                          |   2          2        | 
                          \ dx         dr         / 

The portion of the time used to construct the matrices and vectors is insignificant compared to the portion used for the integration. To see that this is true, in your initialization section include the line

unprotect(int):  int:= 1:

And make the last line of the worksheet 

time() - timer1;

Then execute the entire worksheet using the !!! on the toolbar. Doing this, I get about 1 second. 

Why do you set Digits:= 16? I can see no point in using a value greater than 15 (the cut-off for hardware-float computation) for this computation.

If the matrix and vector construction time were an issue, a small time reduction could be made by replacing 

Matrix(M+1, M+1, (j0,j)-> ...)

with

rtable((M+1)$2, (j0,j)-> ..., datatype= hfloat, subtype= Matrix)

and replacing

Vector(N, i-> ...), et al.,

with

rtable(N, i-> ..., datatype= hfloat, subtype= Vector[column])