An efficient way to do this is to use a little-known Maple data structure called a heap. A heap is like a stack or a queue except that the total ordering of the elements is determined by a user-supplied function rather than by arrival time.

SeqA:= proc(n::posint)
uses NT= NumberTheory;
local
    k, p:= 1, #current and initial value 
    A:= Array(1..1, [p]), #storage for sequence
    #structure that tracks "least entry not previously used this way": 
    H:= heap[new]((i,j)-> `if`(A[i]=A[j], i>j, A[i]>A[j])),
    N:= table(sparse) #non-novel entries
;
    for k to n-1 do
        A(k+1):= `if`(N[p]=0, NT:-tau(p), p + A[heap[extract](H)]);
        N[p]:= (); p:= A[k+1]; heap[insert](k, H)
    od;
    seq(A)
end proc
:    
SeqA(99);
1, 1, 2, 2, 3, 2, 4, 3, 5, 2, 4, 6, 4, 7, 2, 5, 7, 11, 2, 6, 8, 
  4, 8, 12, 6, 11, 16, 5, 11, 16, 22, 4, 10, 4, 8, 12, 19, 2, 9, 
  3, 5, 8, 13, 2, 10, 12, 20, 6, 14, 4, 10, 14, 22, 31, 2, 12, 
  14, 24, 8, 18, 6, 14, 20, 31, 42, 8, 19, 27, 4, 16, 20, 32, 6, 
  18, 24, 36, 9, 22, 31, 45, 6, 20, 26, 4, 18, 22, 36, 52, 6, 22, 
  28, 6, 22, 28, 46, 4, 22, 26, 44

Another way to do it is to set A4(0.):= 1. before doing the seq command. This is called setting a remember table value.

Here is a procedure to convert a piecewise expression into a Fortran-compatible procedure that can be differentiated with D. This uses a feature new to Maple 2021, so if you're not using Maple 2021, let me know.

(*---------------------------------------------------------------------
This procedure converts a piecewise *expression* into a procedure that
    1. can be differentiated or partial differentiated with D,
    2. can be converted to Fortran and its derivatives can be
    converted to Fortran.

The optional 2nd argument is a list of the names (optionally with type
specifiers) that will be used as
the parameters of the procedure, and hence the potential variables of
differentiation. Its default value is the nonconstant names that 
appear in the piecewise conditions with hfloat as the type.

This procedure uses a feature new to Maple 2021.
----------------------------------------------------------------------*)
`convert/pwproc`:= proc(
    P::specfunc(piecewise),
    V::list({name, name::type}):= 
        [indets(indets(P, boolean), And(name, Not(constant)))[]]::~hfloat
)
option `Author: Carl Love <carl.j.love@gmail.com> 2021-May-9`;
    subsop(
        8= hfloat, #procedure return type
        unapply(
            subsindets(
                convert(
                    if op(0,P)::symbol or nops(P)::odd then P
                    else op(0,P)(op(P), op([0,1],P))
                    fi,
                    ifelse, ':-recurse'
                ),
                specfunc(ifelse), convert, `if`
            ),
            V
        )
    )
end proc
:  
#The piecewise expression from your Question:
P:= piecewise(
    x <= 0, x^2+x, 
    x < 3*Pi, sin(x), 
    x^2 - 6*x*Pi + 9*Pi^2 - x + 3*Pi
):
PP:= convert(P, pwproc);
 PP := proc (x::hfloat)::hfloat; options operator, arrow; if x 
    <= 0 then x^2+x else if x < 3*Pi then sin(x) else 
    x^2-6*x*Pi+9*Pi^2-x+3*Pi end if end if end proc

dPP:= D[1](PP); #or simply D(PP)
dPP := proc (x::hfloat) options operator, arrow; if x <= 0 then 
   2*x+1 else if x < 3*Pi then cos(x) else -6*Pi+2*x-1 end if 
   end if end proc

Digits:= 18: #for correct value of Pi as a Fortran "double"
CodeGeneration:-Fortran(dPP);
Warning, procedure/module options ignored
      doubleprecision function dPP (x)
        doubleprecision x
        if (x .le. 0.0D0) then
          dPP = 0.2D1 * x + 0.1D1
          return
        else
          if (x .lt. 0.3D1 * 0.314159265358979324D1) then
            dPP = cos(x)
            return
          else
            dPP = -0.6D1 * 0.314159265358979324D1 + 0.2D1 * x - 0.1D1
            return
          end if
        end if
      end

Here are 5 methods for your map_with_index:

map_with_index1:= (f,a)-> (L-> f~(L, [$1..nops(L)]))(convert(a, list)):
map_with_index2:= (f,a)-> [for local i,x in a do f(x,i) od]:
map_with_index3:= (f,a)-> local i:= 0; map(x-> f(x, ++i), a):
map_with_index4:= (f,a)-> 
    (L-> zip(f, L, [$1..nops(L)]))(convert(a, list))
:
map_with_index5:= (f,a)-> local i; [seq](f(a[i], i), i= 1..nops(a)):
 
L:= ["a", "b", "c"]:
(print@~map_with_index||(1..5))(`[]`, L);
                 [["a", 1], ["b", 2], ["c", 3]]
                 [["a", 1], ["b", 2], ["c", 3]]
                 [["a", 1], ["b", 2], ["c", 3]]
                 [["a", 1], ["b", 2], ["c", 3]]
                 [["a", 1], ["b", 2], ["c", 3]]

What you call flat_map is not worth writing or even giving a name to in Maple. Instead, in the calling of the 5 procedures above, just change `[]` to1@@0.

 
Here are 10 ways to extract a sublist based on a predicate:

(
    remove(`>`, L, "b"),
    remove(x-> x > "b", L),
    remove[2](`<`, "b", L).
    select(`<=`, L, "b"),
    select[2](`>=`, "b", L),
    map(x-> `if`(x > "b", [][], x), L),
    (x-> `if`(x > "b", [][], x))~(L),
    map(x-> if x > "b" then else x fi, L),
    [seq](`if`(x > "b", [][], x), x= L),
    [for x in L do if x > "b" then else x fi od]
);

 ["a", "b"], ["a", "b"], ["a", "b"] . ["a", "b"], ["a", "b"], 
   ["a", "b"], ["a", "b"], ["a", "b"], ["a", "b"], ["a", "b"]

It'll be useful if we can assume that all variables represent positive numbers. If that assumption is not valid, let me know, and I think that some modification of the techique below will still work.

Outline of my steps coded below:

1. Replace all variables with the squares of new variables. The new variables have names beginning z. This step removes all radicals from your equations. They are now equations of rational functions.
2. Subtract one side from the other in each equation. They are now just rational functions implicitly equated to 0.
3. Extract the numerators. Now we have a system of polynomials implicitly equated to 0 (which I called sysP).
4. Use the eliminate command to express everything in terms of the two variables that you mentioned.
5. Replace the new variables with the square roots of the original variables.

sys:= sys1 union sys2:
v:= [indets(sys, name)[]]; 
vz:= subsop~(0= z, v);
sysP:= (numer~@simplify)(subs(v=~ vz^~2, (lhs-rhs)~(sys)), assume= vz::~positive);
#Let's see how complicated those polynomials are.
[degree, length]~([sysP[]]);
           [[2, 101], [2, 109], [26, 418], [26, 414]]

#They're much simpler than I expected, and much simpler than the original
#system! We have 4 short polynomials: 2 of degree 2 and 2 of degree 26.

Sols:= map(
    s-> (lhs^2=rhs^2)~(s[1]) union s[2] =~ 0,
    subs(vz=~ v^~(1/2), [eliminate](sysP, {z[4,1], z[4,2]}))
):
length~(Sols);
                  [35, 122754, 205504, 205504]

We get 4 solutions (and it's done with amazing speed---6 seconds in Maple 2021). The first solution is trivial, and also irrelevant under our positivity assumption. The remaining solutions are extremely long, but they do satisfy exactly what you asked for: expression of everything in terms of y[2,1] and y[2,2]. And, amazingly, this has been done without the use of RootOf.

As of Maple 2019 (or perhaps 2018), embedded for loops, such as you show in your Question, are allowed (in 1D input only!). The syntax is:

n:= 7:
a:= [$1..7]:
P:= piecewise(for i to n do x < a[i], y[i](x) od);

This avoids the op([...]) syntax required by seq.

Yes, you are correct that your piecewise expression can be converted to an equivalent form that doesn't use piecewise and this makes some improvement to the time of the integration. After defining the piecewise via f:= piecewise(...), do

f:= convert(f, Heaviside):

Other recommendations:

Change Digits:= 30 to Digits:= 15, which is the largest value for which hardware-float computations will be used, which are much faster. As far as I can see, this doesn't change the final result, but my observations of this aspect have been superficial thus far.

Remove Digits:= 4. In my opinion, that's too low to ever be used as the global value of Digits. There are other ways to control the precision of a numeric integration, which are discussed in the next paragraph.

Inside the int command, after keyword numeric, add the options epsilon= 0.5e-4, digits= 7. The returned result will be guaranteed to 4 significant digits. There is some chance that the integral will be returned unevaluated, which means that it wasn't able to control the accuracy sufficiently. If this happens, let me know.

If you make these changes, then I think that the majority of the computation time is consumed by solve, not by int. If this amount of time is still too much for you, let me know.

The sequence of iterates in Newton's method can be extremely sensitive to the number of digits of precision at which the calculations are performed. This is totally normal.

You may already know this, but that is not conveyed in your Question or worksheet. Since the second derivative of exp(x) never changes sign (in other words, there are no inflection points), there can be at most two intersections between the curve and any straight line, and if an intersection is also a point of tangency, then it's necessarily the only intersection. The cases of 0 intersections and 1 intersection are relatively easy to classify, as has already been done. Thus, all remaining cases are 2 intersections.

Do you mean this?:

data:= convert(
    Statistics:-Sample(Uniform(-2, 2), [100, 2]), 
    list, nested
):
plot(
    `[]`~(data), 
    color= COLOR~(
        HUE,
        .5 -~ argument~(-(Complex@op)~(data)^~2)/~(2*Pi)
    ),
    style= point, symbolsize= 20, symbol= solidbox
);

Maple has two packages for parallelization: Threads and Grid:

• With Threads, the procedures running on different cores share all memory except for variables that you explicitly tell them not to share. 
• With Grid, the procedures in different cores run as distinct "processes", each with its own memory, and they can only communicate with each other (if any communication is needed at all) via variables or messages that you explicitly send between them.

The vast majority of Maple's high-level symbolic library commands make extensive use of global variables, which makes them unsuitable for the shared-memory environment of Threads. I suspect that Query is in this category.

Try using the command Grid:-Map for this.

I'm guessing that you might want some tutorial help with doing this problem "by hand".

The first step is to "complete the square" in the denominator, thus obtaining 1/((s+4)^2 + 4). We recognize this as the composition f@g of the two functions g:= s-> s+4 and f:= s-> 1/(s^2+2^2). The inner function g is a "shift": just adding a constant a to s. That means the inverse transform will contain exp(-a*t) as a factor (not as a function composition). The remaining part is of the form 1/(s^2+b^2). From a table of transforms (even the most-basic tables will likely have this one), the inverse transform of that is sin(b*t)/b. This is the other factor, making the answer exp(-4*t)*sin(2*t)/2.

A second approach to the problem is much more straightforward---it doesn't require much recognizing of things as being of certain special forms---but it does delve into complex numbers. We factor the denominator over the complex numbers to obtain (s+(4+2*I))*(s+(4-2*I)). Compute the partial fraction decomposition of the rational function from that factorization to obtain 

I/4*(1/(s + (4+2*I)) - 1/(s + (4-2*I))

Each fraction is a "shift", so the inverse transform is

I/4*(exp(-(4+2*I)*t) - exp(-(4-2*I)*t))

If you convert that to real form (which can be done with Maple command evalc) you get the same thing as by the previous method.

Your transcription of your code is so chock full of typos that the only reason that I can understand your Question is because I worked on an earlier Question of yours to which this one is closely related. I suspect that very few, perhaps 0, other readers can understand it at all, which may be why it's taken a long time for you to get an Answer.

In addition to the typos, the expressions -Pi/2 + 10 and Pi/2 - 10 are surely wrong. By 10, did you perhaps mean 10 degrees? In that case, you need to change 10 to Pi/18.

Now let's address your undefined issue: As is well known and explicitly taught---even in first-semester calculus---the derivatives of piecewise functions are often undefined at points where one piece changes to the next. So the fact that you're getting undefined in your piecewises is not an error or a bug (and I'm not saying that you were claiming that it was an error or a bug). But I understand that some further symbolic processing that you're doing with these expressions is having problems because of the undefineds. That's totally believable. The following example serves to illustrate where these undefineds come from, and also shows a method to avoid them. Please execute the code:

f1:= abs(x):
f2:= convert(f1, piecewise, x);
diff(f2, x);
#Note the correct presence of undefined in previous result.

f3:= convert(f2, Heaviside);
convert(diff(f3, x), piecewise, x);
#Same result as before

Heaviside(0):= 0:
convert(diff(f3, x), piecewise, x);
#Now there's no undefined.

There are many subtleties involved in working symbolically with piecewise expressions, especially when there is more than one independent variable (it looks like you have both x and y), so I'm not claiming that the above will be a complete solution to the undefined issue. See if you can make *some* progress with it, and please report back.

I agree that the result 1 doesn't make sense and should be corrected. But compare your result with this:

convert(abs(1,x), Heaviside) assuming x::real;
                      -1 + 2 Heaviside(x)

This small example, as well as Tom's help page excerpt, should answer your Question about $:

f1:= proc(x, y) x^2+y^2 end proc:
f1(a, b, c);
                             2    2
                            a  + b 

f2:= proc(x, y, $) x^2+y^2 end proc:
f2(a, b, c);
Error, invalid input: too many and/or wrong type of
arguments passed to f2; first unused argument is c

The x and y in those examples are called parameters, not arguments. The a, b, and c are the arguments.