@Carl Love I should've included some explanation of my plot command. So, here it is: I created two separate plots and overlaid them with plots:-display. The first of the inner plot commands plots the shaded region; the second plots the curve. Any number of plots in the same dimension (created with any commands, not just plot) can be overlaid with plots:-display.

A single plot command can also handle multiple plots, thus avoiding the need for display. However, if the x-ranges are different, the way that I showed is easiest. But finer control---allowing for both multiple x-ranges and vertical lines within a single plot command---can be achieved by specifying the curves parametrically (i.e., specifying their x- and y-components as separate functions):

plot(
   [
      [x, f, x= 0..1], #shaded (filled) region
      [x, f, x= -1..1], #full curve
      [1, y, y= 0..eval(f, x=1)] #vertical line segment
   ], 
   filled= [true, false, false], thickness= 3,
   labels= [x, ``], labelfont= [TIMES, BOLD, 16],
   legend= [A, f, x=1]
);

If you want to be very fancy, labelling text can be put directly in the plot:

plots:-display(
   plot(
      [
         [x, f, x= 0..1], #shaded (filled) region
         [x, f, x= -1..1], #full curve
         [1, y, y= 0..eval(f, x=1)] #vertical line segment
      ], 
      filled= [true, false, false], thickness= 3,
      labels= [x, ``], labelfont= [TIMES, BOLD, 16]
   ),
   plots:-textplot(
      [
         [0.75, 0.5, typeset(A), font= [HELVETICA, BOLD, 24]],
         [
            0.5, eval(f, x= 0.5), typeset('f'=f),
            #make baseline of text parallel to f's tangent line: 
            rotation= arctan(eval(diff(f,x), x= 0.5)), align= above
         ],
         [1, 1, typeset(x=1), rotation= Pi/2, align= right]
      ],
      font= [TIMES, BOLD, 16]
   )
);

@Oliveira There's a bug in Acer's example, which is totally excusable since he was posting from his phone. The &-operators have a higher precedence (i.e., their position in the order of operations) than all other arithmetic operators. Thus, any arguments containing arithmetic operators need parentheses:

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

Acer's version cannot handle multi-argument commands. I'll leave it to him to modify that, since his style may avoid parameterless procedures, such as I used.

In what computer/mark-up language is this expression (copied directly from your Question) written?:

∫(-625 R^2 ro (-2 cp3 x1 lambdaopt+sin((pi t)/10) cp2+(16 cp2)/5) (-cos((pi t)/10)+cos(2 pi)+((-t/5+4) x1+(8 t)/25-32/5) pi) lambdaopt pi b11 (e)^(-((16+5 sin((pi t)/10)) cp1)/(10 lambdaopt x1))+625 (-R^2 k11 (cos((pi t)/10))^3+k11 R^2 (cos(2 pi)-((t-20) (x1-8/5) pi)/5) (cos((pi t)/10))^2+((32 sin((pi t)/10) R^2 k11)/5+(281 R^2 k11)/25+4 lambdaopt^2 (a11 x1+a13 x3)) cos((pi t)/10)-(32 k11 R^2 (cos(2 pi)-((t-20) (x1-8/5) pi)/5) sin((pi t)/10))/5+(281 (x1-8/5) (R^2 k11+(100 lambdaopt^2 (a11 x1+a13 x3))/281) pi t)/125) x1^2)ⅆt

@Oliveira There's a long-standing bug in the kernel regarding eval that forces me to make a small change to the operator to handle your example. The new operator is

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

Then your eval example would be handled by 

(x^2 - 2*x, x=4) &( eval

@tomleslie The problem does explicitly say that the function should be plotted from -1 to 1, not merely that that be the viewing window.

@Samir Khan Ah, yes, I figured something else was coming. Building the aniticipation, eh?

@dharr But Heaviside(x) <> x throughout any real open interval.

Surely you understand the very limited value of this Post if you don't reveal the process of extracting the equations from from the cursive script, which surely preceded the equations. 

@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