@hamideh The are a great many commands in Maple that accept equations or algebraic expressions interchangeably. The algebraic expressions in these cases are treated as if they were equated to 0. If A and B are algebraic expressions and R is a boolean expression, then `if`(R, A, B) is still an algebraic expression; but if E is an equation, then `if`(R, E, B) is not an equation because the primary (or outermost) operator of an equation is `=`.

What definition of continuous are you using for your students? Even if everything is restricted to real numbers, I'd still say that sqrt is continuous: The set of squares of any open interval of the codomain is (relatively) open in the domain [0, infinity). Being undefined for some real arguments isn't quite the same thing as being discontinuous.

@Muhammad Usman And I suppose that you need me to tell you how to make that change? Obviously, an M2 needs to be changed to an M1 somewhere in the code. There are only 3 M2s in the code. How difficult can it be?

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?