@mmcdara Yes, I figured that they were node values of a numeric method. If so, they could probably be easily derived. One surprising and potentially problematic aspect is that all 6 of the negative values have significantly lower magnitude than all 6 of the positive values. The most extreme such ratio is approximately 900 (just eyeballing it), or nearly 3 digits. This can lead to a round-off error problem sometimes called catasthropic cancellation.

@DimaBrt It may be possible to fine-tune Maple's built-in numeric ODE BVP solver to manage memory usage for your large system. Some options that come to mind are

1. method= bvp[trapdefer] (default is traprich)
2. maxmesh= 2^5 (default is 2^7)
3. abserr= 1e-3 (default is 1e-6)

See help page ?dsolve,bvp,numeric.

If this was your homework, and I was grading it, the first thing that I would say is your array b is a list of magic numbers---complicated numbers of mysterious origin that appear in a program without explanation or derivation. At the very least, you need to comment on what they mean, where they came from, and how they were derived. Better yet is to derive them in the program (if that's possible) or extract them from some source document (if it's not).

@vv Regarding the "oddness" of the integral: Nobody has addressed my concern that this integral appears to be missing a change-of-coordinates factor.

@nm 

The command hastype is just a specialization of the much more frequently used command indets, which actually extracts the subexpressions rather the just saying true/false about their presence. So, hastype(e,T) is essentially^[*1] equivalent to evalb(indets(e,T) <> {}). (Note that indets always returns a set, possibly empty.) So, the way to achieve your goal is to use patmatch (or perhaps typematch) to parse the element(s) of the set returned by indets. This is a much easier way to use patmatch than what you were trying to do because you can still completely ignore the surrounding expression(s) in which your target(s) are embedded. Note that this technique will work just as well for multiple targets or multiple expressions.

I'll write some code for your example in several hours. The possible presence of a numeric coefficient a0 (other than 1) is a complicating "factor" that nonetheless can be dealt with; however, if you really actually don't care about a0, let me know, because by ignoring it, the complexity of the code can be cut in half (or so). 

[*1] I say "essentially" because hastype is built-in to the kernel, so the actual mechanism by which it works may not be exactly as I described, whereas the final results are.

@nm When using types (as with hastype, indets, subsindets, etc.) you only need to encode a generalized form of the target expression; you can completely ignore encoding anything about the surrounding expression in which the target is embedded. That's a major reason why hastype is much easier to use than patmatch.

There's no limit to the complexity of a type: Anything that can be encoded as a boolean check can be used^[*1]. Procedures can be used^[*1]. Types can be named and stored with TypeTools:-AddType. Once named, types can be used recursively or used to build more-complicated types. 

[*1] Just make sure that your boolean checks and procedures will return true or false for any single-expression input, no matter how weird. They should never return error, FAIL, or not fully evaluated.

@mmcdara 

Wolfram MathWorld^[*1] defines the "sample central moment" via the same simple (biased) formula used by Maple's CentralMoment. (And note that the mean used in the formula is the sample mean---so, clearly, it has lost a degree of freedom, and is thus biased.) At the bottom of the page, the bias correction factors for moments 2-5 are given.

I can see a good argument that Variance should offer an option to specify sample or population. At the very least, the formula being used must be mentioned on the help page (likewise for StandardDeviation). But I think that the central moment is a more "primitive" or "raw" concept, so CentralMoment should remain as is.

[*1] Weisstein, Eric W. "Sample Central Moment." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/SampleCentralMoment.html

@minhthien2016 I don't mean any disrespect toward Mariusz's Answer; indeed, I gave it a vote up. But isn't VV's technique far superior? With your revised problem, I think that simplify gets bogged down on the hugely complicated algebraic numbers produced by solve. I didn't wait for it to finish.

@nm I didn't think that has would actually be enough for your overall purpose. But I do think that hastype is enough:

hastype(expr, identical(x*y)^And(fraction, satisfies(e-> 1/6 <= e and e <= 1/3)));

I would strongly persuade anyone toward hastype and away from patmatch in almost all cases. It's simply much easier to construct the patterns (aka structured types).

@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.