The following little procedure is a direct one-for-one replacement for the deprecated linalg:-blockmatrix:

BlockMatrix:= proc(m::posint, n::posint, L::list({Matrix, Vector}))
local k;
   Matrix([seq(L[(k-1)*m+1..k*m], k= 1..n)])
end proc:

@dharr The information that can be obtained that way is incomplete. One would hope that Statistics:-Distribution(X) would return Normal(0,1), but it instead returns 

Error, (in Statistics:-Distribution) invalid input: too many and/or wrong type of arguments passed to NewDistribution; first unused argument is _R
 

@ecterrab Note that a distinction needs to be made between the index being NULL and there being no index at all. In the call f[](a), the index exists, but is NULL. The author of f may wish to distinguish this from the call f(a), which has no index at all.

While this distinction may seem pedantic, it is actually used in piecewise. Note the difference between

f1:= piecewise(x > 0, 1);

and

f2:= piecewise[](x > 0, 1);

AFAIK, this useful feature of piecewise is not documented.

D is an example of a standard library procedure such that D(...) and D[](...) are both meaningful and those meanings differ.

@Muhammad Usman Okay, what you want now truly is a block matrix, whereas the prior thing wasn't. It can be created like this:

restart:
(M1,M2):= (3,3):  nu:= 1:
F:= Matrix(
   convert(
      Matrix(
         (M1,M2), 
         (r,s)-> Matrix(
            [[`if`(r > s and (r+s)::odd, 4*(s-1+nu), 0)$M2]],
            shape= diagonal, scan= band[0]
         )   
      ),
      list, nested
    )
);  

The purpose of my outer two commands---

Matrix(convert(..., list, nested))

---is to convert it from an explicit block matrix to a normal (flat) matrix. These commands effectively replace the deprecated linalg:-blockmatrix.

@Kitonum In Maple 2019, a neutral operator can be used in functional form without quotes. Prior versions require quotes.

&x(1,2)  #works in Maple 2019, but not earlier
`&x`(1,2)  #works in any 
1 &x 2  #works in any

@DJJerome1976 The correct word is co-domain rather than domain, but I understand what you mean. And yes, it can be restricted to a real co-domain like this:

f:= piecewise(x >= 0, sqrt(x), undefined);
iscont(f, x= -infinity..infinity); #Domain is real numbers.
iscont(f, x= 1..10); #Domain is bounded real interval.

@Muhammad Usman What do you mean by "any value of M"? The code I already gave will work for M being any nonnegative integer. Isn't that good enough?

@Mariusz Iwaniuk I suspect that the above code will only work in 1-D input (aka Maple Input).

@erik10 You're welcome. Both of the above can be generalized to Arrays of any number of dimensions or to switch the roles of rows and columns.

I may be missing some crucial detail here, but doesn't the symbolic solution provided in Rouben's Answer to your previous Question also completely answer your new Question? That solution is

x[2](t) = (p[2]+1)^2
*p[1]*cosh(t)^2+(-2*(p[1]*p[2]^2-p[3])*(p[2]+1)*p[1]*sinh(t)-p[2]^2*p[1]-2*p[1]
*p[2]-2*p[1]+p[3])*cosh(t)+(p[1]*p[2]^2-p[3])^2*p[1]*sinh(t)^2+(2*p[1]^2*p[2]^2
-2*p[1]*p[3]+p[2]+1)*sinh(t)+p[1]
 

@Rouben Rostamian  By setting Digits:= 15 and using dsolve options abserr= 1e-10 and relerr= 1e-10, I got the accuracy I wanted. I didn't need to fiddle with the events, nor could I figure how to do so. That help page is quite unhelpful, having about 10 times as many options as there are examples. So, combining the higher accuracy with your way of programmatically extracting the values, I am satisfied.

@Rouben Rostamian  I see two drawbacks to this approach:

1. Its default accuracy is much, much worse than the other "stock" methods, but I guess that I could figure out how to adjust that with the relerr option to dsolve(..., numeric).
2. As you've written it, there's no way to programmatically access the returned valued. Do you know a solution for that?

@Carl Love Since the u's only appear with phi applied to them on the right side of the defining equation, if phi is injective over any relevant domain, then without loss of generality (WLOG) we can assume that phi is the identity function. (And if phi isn't injective, then I think that the problem is ill posed.) So, phi is just obfuscation, and the notation can be simplified to

`&*`:= (u__1,u__2)-> ln(exp(u__1)+exp(u__2)) + Q(u__1 - u__2);

@sunit You should read the help page for algsubs (?algsubs). It provides numerous examples of some of the subtleties involved. It is easy to learn to recognize syntactic subunits, and once you can, there's no question that for them you should use eval (which understands mathematical subtleties) or subs (which intentionally ignores them) to make the substitutions. In particular, any name (as in a variable name) is always a syntactic subunit.

There are many commands for making substitutions. Here they are ordered by my approximation of their complexity (and I've probably forgotten a few). I'd generally recommend choosing the simplest one on the list that you know will get the job done.

subs
subs[eval]
eval
eval[recurse]
subsindets[flat]
subsindets
evalindets[flat]
evalindets
simplify(..., {...}) #aka "simplify with side relations" (see Kitonum's Answer)
algsubs
eliminate

If you have a question about factor, you should ask that as a separate Question.

Are you saying that you want to use Maple's (G)UI as an interactive front end for a Linux shell session with the two communicating both way via pipes? Or would it suffice to just use Maple as an editor to compose (and perhaps to analyze) a batch of shell commands that could be processed by Linux as a batch?