The difference is simply that CentralMoment is treating the data as a whole population and Variance is treating it as a sample. To me, it seems useful to be able to make that distinction, so I like that the two commands return different results.

You wrote:

• By definition the variance of a random variable X is its 2nd order central moment. This definition should also apply when these statistics are calculated on a sample of X.

But how is the command supposed to know whether the data that you've given it is a sample or the whole population?

Here's a procedure with an optional second argument to specify the depth of de-nesting. If the second argument is omitted, it does exactly what Acer's code does.

Flatten:= proc(S::set, depth::nonnegint:= infinity)
local s, r:= S;
   if depth=infinity then return {subsindets(S, set, op)} fi;
   to depth do r:= {seq(`if`(s::set, s[], s), s= r)} od;
   r
end proc
:

Examples:

A1 := {1,2,{5,6},3,{13,{{21,27},16}},{8,9},{{11,12},4}};
Flatten(A1); Flatten(A1,1); Flatten(A1,2);
 

Do 

remove~(evalb, S)

The image sqrt(R) is a connected set in the complex plane^[*1]: It's the union of the nonnegative real and the nonnegative imaginary semi-axes. Just because it's not real-valued for all real arguments doesn't automatically make it discontinuous. 

[*1] By the way, I'm not claiming that a connected image necessarily implies continuity; it doesn't. But by considering this image, it's easy to see that sqrt is an embedding of R into C and hence continuous.

I don't see any practical difference between what you ask for and 

A:= Array((0..M)$3, ()-> a[args])

Rather, the only difference that I see is that you call it a "block matrix", whereas I call it a "three-dimensional array".

It seems like it would take a great many steps to convert an ODE from Maple standard form (i.e., the form(s) accepted by dsolve) to the specialized form required by Edgardo's suggestion, DifferentialAlgebra:-Tools:-LeadingDerivative. I look forward to seeing an example of it. On the other hand, it looks like it'll work for PDEs also.

Here's a procedure to do it just for ODEs:

ODEdegree:= (e::{algebraic, `=`(algebraic)})-> 
   try
      (e-> degree(
         e, 
         indets(
            e, 
            typefunc(
               (O-> `if`(O>1, typefunc('`@@`'(identical(D),O)), specfunc(D)))
                  (PDEtools:-difforder(e))
            )
         )
      ))
         (evala(Norm(convert(`if`(e::`=`, (lhs-rhs)(e), e), D))))
   catch:
      FAIL
   end try
;

It works for both equations and expressions, and for D-form and diff-form (or even a mixture). Here it is applied to your example:

ode:= (1+diff(y(x),x)^2)^(3/2) = diff(y(x),x$2);
ODEdegree(ode);
      2

It looks like evalc won't automatically assume that D(x11)(1) and D(y11)(1) are real. It will assume that if diff form is used instead of D. So, this gets you your desired result:

convert(evalc(Re(convert(F1,diff))), D)

If you plan on doing a lot of this, I could add code to evalc (specifically, a procedure named `evalc/D`) that would do the above automatically.

There seems to be a critical-severity bug in fsolve in Maple 2019.1 such that every attempt that I've made at solving this with fsolve (using various formulations of the function whose root is wanted) has led to a kernel crash. This is particularly disturbing in this case because the function is differentiable, strictly increasing, and has a unique real root.

Another root-finding command--an alternative to fsolve--is RootFinding:-NextZero. Here it is applied to your problem. I've done this in a way that can be applied to most arclength problems:

restart:

#These first 2 steps are simply to reconstruct your original function
#from its derivative:
C:= 5.557990765:
f:= int(<-C,-7.3,-5.6,7.3>.<[sin, cos, sinh, cosh](C*x)>, x);

x0:= 0: #Start of arclength measurement.
L:= 0.5: #Desired arclength.

F:= subs(  #Construct numeric procedure whose root we want.
   [_ArcLInt= VectorCalculus:-ArcLength(<x, f>, x0..X, 'inert'), _ArcL= L],
   proc(X) local x; evalf(_ArcLInt) - _ArcL end proc
):    
x1:= RootFinding:-NextZero(F, x0);
                       x1 := 0.1550088023
plot(f, x= x0..x1);

An easy way to get a close upper bound on the desired x-value is to use the straight-line distance between two points on the curve:

x_ub:= fsolve((x0-x)^2 + (eval(f, x= x0) - f)^2 = L^2);
                      x_ub := 0.1570169605

This upper bound might be useful if you're using a root finder whose performance is improved by specifying bounds.

An example of an fsolve that leads to a kernel crash is

fsolve([F], [x0 .. L])

using procedure F and x0 and L from above.

To permute the rows of Matrix M so that they're sorted (ascending) by column C, do

M[sort(M[..,C], output= permutation), ..]

To change that to descending, make `>` the second argument to sort.

To count the entries of column C that are larger than V, do

add(`if`(x > V, 1, 0), x= M[..,C])

You can read the through the code using showstat(fsolve), showstat(fsolve::sysnewton), showstat(fsolve::linsolve), etc.

One eye-opener for me is that linsolve uses a directly coded Gaussian elimination rather than using the extensively optimized LinearAlgebra:-LinearSolve. This is done for every iteration of Newton's method.

It only evaluates when you press RETURN. If instead you end lines with SHIFT-RETURN, you'll get a new line without any evaluation happening. You can enter any  number of lines this way. Each group of such lines is called an execution group.

You need to use algsubs instead of subs:

eqn:= N*x*y - x*y + f(x,y) = 0;
algsubs(x*y= 1, eqn);

Both subs and eval (and even subsindets and evalindets) can only substitute for subexpressions that occur as syntactic subunits in the main expression, but algsubs can find algebraic subunits. It would require a lengthy digression into the minutiae of computer algebra for me to precisely define syntactic subunit, but if you'd like, I'll try later. I think an example will suffice as a definition for now: Let's consider the present expression, but written in functional form (I'd also call it a tree form):

`=`(`+`(`*`(N, x, y), `*`(-1, x, y), f(x,y)), 0);

The above is not just a pedantic presentation form, but is actually executable Maple code that produces the same eqn as the original. You will note that our target expression `*`(x*y) doesn't appear in the above functional form, which is why subs doesn't see it.

The above functional form is very close to (but not identical to) Maple's internal representation of the expression. In reality, the numeric coefficients are handled a little bit differently, which is why subs does find the x*y in `*`(-1, x, y).

For the purpose of checking associativity, we can ignore the phi on the left side of the equation. Then do this:

`&*`:= (u__1,u__2)-> ln(exp(phi(u__1))+exp(phi(u__2))) + Q(phi(u__1) - phi(u__2));
AssocLeft:= (u__1 &* u__2) &* u__3;
AssocRight:= u__1 &* (u__2 &* u__3);
CriticalEqn:= AssocLeft - AssocRight = 0;

So, to check any particular Q, we substitute it into the CriticalEqn, along with some arbitrary phi and u's. Here, I let Q and phi be the identity function:

eval(CriticalEqn, [Q= (x->x), phi= (x->x), u__1= 1, u__2= 0, u__3= -1]);
evalf(%);

The left side isn't even close to 0, which shows immediately that it's not associative for this Q.

Regardless of whether the equations are linear or the solutions are exact, you can achieve what you asked for like this:

eval(Vector[row](nops(SOL[0]), symbol= tau), SOL[0]);

Here is code to exactly duplicate in Maple 2019 the functionality of rept that you show in your worksheet:

rept:= (f, R::(name= range(realcons)))-> 
   ((t,a,r)-> eval(f, t= t - r*floor((t-a)/r)))((lhs, (lhs, rhs-lhs)@~rhs)(R))
: 

Your example:

F:= rept(piecewise(t<1, t/2+1/2, t>=1, -t+2), t= -1..2);
F := piecewise(t - 3*floor(t/3 + 1/3) < 1, t/2 - (3*floor(t/3 + 1/3))/2 + 1/2, 1 <= t - 3*floor(t/3 + 1/3), -t + 3*floor(t/3 + 1/3) + 2)

Note that my code gives an exact symbolic representation (in algebraic form) of the periodic extension, not merely a procedure or a function only suitable for plotting.