The problem doesn't really ask you to plot the line x=1. If it did, there are many ways to do it, one of which (implicitplot) is shown by DJJerome, and another of which (parametric plot) is shown in my Reply below. It being a vertical line complicates the syntax a tiny bit. Here's a way that just plots the function over the stated interval with the area A highlighted:

f:= x*exp(x^3):
plots:-display(plot(f, x= 0..1, filled), plot(f, x= -1..1));

[I was writing this while Acer was writing his similar Answer.]

It's trivial to create a custom binary (i.e., two-argument) operator^[*1] to do what you asked for:

`&(`:= ()-> args[-1](args[1], args[2..-2]):

You might as well put the above in your initialization file (or in a custom library if you have one). Its usage is thus:

w2:= (sin(x)^2 + cos(x)^2) &( simplify;

or, if you wish to assign the left-side value also:

w2:=  (w1:= sin(x)^2 + cos(x)^2) &( simplify; 

Multiple arguments can be accomodated by putting them in parentheses on the left:

(x^2*y^3, y) &( degree;

Even when the left side is a single argument, it needs to be parenthesized in most cases.

Chains (compositions) of operations can be accomodated:

(sin(x) + cos(x))^2 &( expand &( simplify;

The name &( is my choice. You could make it almost anything that begins with &. See ?neutral for the complete rules. Many, including &(, won't work in 2-D Input, which I never use. But &a, etc., will work.

[*1] Maple calls these neutral operators. They can be used with or without definitions, but, to my mind, once they have definitions they're no longer neutral. So, I called them "custom binary operators".

Use

eval(S1, t= 0.3)

or

S1P:= unapply(S1, t):
S1P(0.3)

The unapply turns it into a procedure and only needs to be done once. After that, the more-convenient S1P(...) syntax can be used for the rest of the session instead of eval(S1, t= ...). Also, the unapply localizes the t so that things won't get messed up if it's inadvertently assigned a value at some point.

To check for a specific function, use type specfunc(...)---in this case specfunc(`odetest/PIECEWISE`). The type function is extremely common, and it's certainly not robust to check for it. For example, sin(x) is a function.

I'm not familiar with DTM, but the ms in both sums in the for loop seem suspicious to me, including the upper limits. I think that all 6 should be changed to k.

Getting Maple to do those computations is trivial, and there are many ways to do it. I show one way below. But there'll be great variation in the efficiency of those ways depending on how the input is supplied, and how your raw computation sequence is converted to a Maple procedure.

Getting Maple to write out a CSV file is also trivial, and it's done with the command

ExportMatrix(file name, matrix name, 'target'= 'csv') 

(see help page ?ImportMatrix).

But I'm confused as to why you "have this file", which contains both one instance of the input data and the computation sequence.

1. Would it be necessary for Maple to read and interpret the computation sequence from this file?
2. Or could you just manually place the computation sequence in a Maple worksheet?
3. Is it necessary that the input data come from such a file?
4. Or could it come from a CSV file?

Suppose that your answers to 2 & 4 are "yes". Then you could do this:

restart:
MyFunc:= evalf@unapply(
   eval[recurse](
      <T | U>, #Output a 2-element row vector.
      #One of many ways to specify computation sequence:
      [T= E^2+G-sin(F), U= F+E-B, E= A*B, F= cos(A)+exp(B), G= 265/C]
   ), 
   [A,B,C] #Input parameters
):
Process:= proc(In::Matrix)
local n:= upperbound(In)[1], Out:= Matrix(n,2), k;
   for k to n do Out[k,..]:= MyFunc(seq(In[k])) od; #Process input row by row.
   Out
end proc
:
Data:= ImportMatrix( some file name , 'source'= 'csv');
ExportMatrix( some other file name , <Data | Process(Data)>, 'target'= 'csv');

If the answers to 1 or 3 are "yes", it's still possible, but I'll need details.

Is there anything wrong with 

has(expr, sqrt(x*y))

After 

infinies:= Infinitesimals(PDE)

do something like

infinies:= eval(infinies, [_F4= 0, _F3= 0, _F2= 1, _F1= (R-> a*R)])

Then, as you said, the Invariants will be much simpler; indeed, they'll fit on one line and contain no unevaluated integrals.

Note 1: Anywhere in Maple, an explicit numeric constant (such as 0 or 1 above) can be used as if it were a constant function.

Note 2: The underscore characters are essential parts of the names _F1, etc., just as if it were an alphabetic letter.

Warning: The _Fs will not necessarily be numbered the same every time that you run this. I advise that you visually check it each time.

There are three problems that I see immediately, before even running your code:

1. In your for loop, the four expressions on the left sides of the assignment operators each have an extra space after the function name. In others words, you need to change S (k+1) to S(k+1), etc.
2. On the right sides of the assignment operators, you need to change all square brackets [ ] to parentheses ( ). Square brackets are never used for algebraic groupings in Maple (nor in most other computer languages).
3. N needs to have a value.

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