It's a stochastic matrix, therefore its eigenvalue of largest magnitude is 1 and all other eigenvalues have magnitude <= 1. See this Wikipedia article: Stochastic matrix

It can be done like this:

sys:= {bcs||(1..2), ode||(3..4)}: 
SFvsNU:= Nu->
    eval(
        diff(f(eta), eta$2),
        dsolve(eval(sys, [N= Nu, a= 1, S= 1]), numeric)(0)
    )
:   
plot(SFvsNU, 1..10);

Both jobs can be done by the same command:

collect(numer(B(4)), x, factor);
collect(F(5), x, factor); #F(5) could also be numer(F(5));

If you want to keep the p factors in each term grouped together, then

collect(numer(B(4)), x, ``@factor);

In the future, please attach a worksheet (using the green uparrow) or include plaintext code of your expressions such that we can cut and paste it into Maple.

You can move mod to the end, like this:

solve(5*x=7, x) mod 16;

That'll give you simply 11.

You can modify mod itself so that it always returns an expression for "all solutions" like this:

restart:
modp__orig:= eval(modp):
unprotect(modp):
modp:= proc(e,p) 
local r, Z:= `tools/genglobal`(_Z);
    assume(Z, integer);
    r:= modp__orig(e,p);
    try r+p*Z catch: r end try
end proc:
protect(modp, modp__orig):

solve(5*x=7, x) mod 16;
                          11 + 16 _Z0

If you're using mods instead of the default modp, then replace modp with mods in all the above.

Like this:

deg:= n-> cat("", n, "°"):
plots:-polarplot(
    t, t= 0..2*Pi, 
    axis[angular]= [  #or axis[2]
        tickmarks= 
            [seq](
                k*Pi/16= `if`(irem(k,4)=0, typeset(deg(k*45/4)), ""), 
                k= 0..31
            ), 
        gridlines= [32, majorlines= 4]
    ]
);

I don't think that there's any good reason why the angle axis is axis[2] rather than axis[1]. That's just the way it is. [Edit: Or you can use axis[angular] instead of axis[2], but that's limited to when axiscoordinates= polar. If using axis[2], then the options work for any 2D plot.]

[Edit: Of course, Kitonum's angularunit= degrees is simpler, but if you want the degree symbol (superscript small o), then I think that you still need the tickmarks option. And since the tickmarks option changes the gridlines, to restore the normal gridlines, I think that you need the gridlines option also.]

I will answer your questions about the scope rules, but there's two preliminary issues to clear up first:

Packages vs modules: Packages and modules are not exactly the same thing. The rules below apply to all modules, regardless of whether they are also packages. The vast majority of modules are not packages, so it would be best if you didn't confuse the issue by using the word "package" for this.

local vs global: Here I am contrasting only the keywords local and global rather than contrasting the concepts of "local variable" and "global variable". It's natural to think that these keywords do parallel things, each "creating" a type of variable. However, that isn't true. In my opinion, global should not be used. The alternative is shown below. I'm not saying that global variables shouldn't be used, just that the keyword global shouldn't be used.

Scope rules:

You asked:

• Is the local variable XM accessible to the procedures used in SUBPACKAGE1?

Your tentative answer for this one was wrong. Yes, it is accessible simply as XM, no matter how deeply nested the procedures and submodules are.

• Is the local variable X1 accessible to the procedure used in MAINPACKAGE ?

You are correct that it isn't (ordinarily) accessible. (For debugging purposes it's possible to make it accessible by setting an option that I won't mention here.)

• Is the global variable YM accessible to the procedures used in SUBPACKAGE1? 
  Is the global variable Y1 accessible to the procedure used in MAINPACKAGE ? 

Global variables are accessible anywhere, regardless of whether they're declared with global. If there's a need to distinguish between a local x and a global x, or if an assignment is being made to global x, use the syntax :-x.

• Is the export variable ZM accessible to the procedures used in SUBPACKAGE1? yes
  Is the export variable Z1 accessible to the procedure used in MAINPACKAGE ? yes with the syntax SUBPACKAGE:-Z1

Your tentative answers to those two are correct; however, for the second question, it makes me cringe that you call it a "package" instead of a module. It would be best if you simply forgot about packages, at least until you've mastered modules.

Your matrix P isn't necessarily stochastic, i.e., its rows don't necessarily sum to 1. You can make it stochastic by redefining some variables, such as

pi[1,2]:= 1-pi[1,1]:
pi[2,2]:= 1-pi[2,1]:

That'll take care of the first two rows. Rows 3 & 4 have many more variables. Did you really intend to use both rho and p as variables?

If you make the matrix stochastic in this manner, your procedure will take about 30 seconds.

Another issue: This is apparently not a problem in your procedure, but I'd recommend against using the use command in a procedure. Instead, use a uses clause, as in

steadyStateVector:= proc(P::Matrix(square))
uses LA= LinearAlgebra;
local n:= upperbound(P)[1];
    LA:-LeastSquares(<P - LA:-IdentityMatrix(n) | <1$n>>^+, <0$n, 1>)^+
end proc
:

or

steadyStateVector:= proc(P::Matrix(square))
uses LinearAlgebra;
local n:= upperbound(P)[1];
    LeastSquares(<P - IdentityMatrix(n) | <1$n>>^%T, <0$n, 1>)^%T
end proc
:

You say that you tried the unassign command and that it didn't work. I suspect that you tried

unassign(phi);

However, (and unfortunately) the correct syntax is

unassign('phi');

A more commonly used alternative to unassign has already been mentioned:

phi:= 'phi';

(Indeed, given that the unassign command requires the quotes, I don't even see the point for its existence.)

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.