@sand15 

The code that you claim (in the Reply immediately above) was proposed by Vv in his Answer---

restart: 
with(Statistics): 
N :=3;
X := RandomVariable(Binomial(N, 1/2)): 
F:=CDF(X, s):

---is not what he actually wrote:

restart:
with(Statistics):
#N :=3;# ... Some integer value >= 2;  
X := RandomVariable(Binomial(N, 1/2)): 
F:=CDF(X, s):
N:=3:

The difference between the two is small typographically, but it is the crucial difference between an obvious bug and its workaround. Please find that difference and play with it until you understand what makes it essentially different. Hint: I already mentioned it in another Reply above.

Don't you think that we check each other's Answers here? I had checked and voted up his Answer within 15 minutes of his posting of it. (That's not to say that I give every correct Answer a vote up---style and/or ingenuity count.)

You see, we're essentially arguing over different pieces of code, which you thought were the same. And who knows?: perhaps you still think they're the same.

Please post your code in plaintext form and/or upload a worksheet. Nobody feels like retyping your code. Most people wouldn't want to Answer a Question if they couldn't test that Answer in Maple.

@acer I also used ToInert to reach exactly the same conclusion as you did.

This is a new low for the 2D Input parser for two reasons:

1. It's impossible to work-around this bug while continuing to use 2D Input;
2. it's impossible to "see" this bug by converting the 2D code to 1D.

I'm not aware of any other bug that has either of those ignoble properties.

Point 2 reveals that the 2D parser operates much more mysteriously than I previously thought. I thought that it first converted the code to 1D (just as if you pasted the code to a 1D field) and then parsed that 1D code the way that 1D code has always been parsed. The 2D parser ought to be required to operate in the manner that I just described! At least that way a savvy user could see what's wrong by converting to 1D, (and an unsavvy user could post it here so that a savvy user could do the same for them). So, this is the last straw for me and 2D Input: Until the red sentence above is "the law", I must say that 2D Input is unsafe to use for any executable code because it is impossible to know with 100% certainty what the actual code even is (which seems contrary to the spirit of mathematics). 

Investigating this, I discovered another violation of "the law": A do .... until ...; statement entered as 2D will parse correctly, but its 1D conversion is do ... until ... end do; which is (unfortunately) a syntax error.

@mmcdara Nobody here doubts that you've found a genuine (and very serious) bug in the CDF of  Binomial---that those commands that you just gave produce a CDF with a serious error!!! So you don't need to send me any file to confirm it.

But these commands, discovered by Vv, produce a correct (although ugly) CDF:

restart:
X:= Statistics:-RandomVariable(Binomial(N,1/2)):
CDF:= Statistics:-CDF(X, s):
N:= 3:
CDF;

This is a workaround. Do you know what that means? It is not a denial of the existence of a bug! It acknowledges the bug's existence and provides the user with an easy and practical way to "work around" the bug until it is fixed.

But I think that Vv is right about you: You're more interested in arguing than in learning practical information about how to use Maple. I noticed that about you in your very first post on MaplePrimes (as sand15).

@mmcdara You totally, totally, missed Vv's point, which is not discont; he just included that as an extra. I suspect that you didn't type his commands into your Maple 2018, because then you'd see that his CDF and its plot are precisely a "piecewise constant [function] with breaks at x=0, x=1, ...x=N".

His point---his workaround---is to first extract the CDF for symbolic N and then instantiate N:= 3, the opposite of the order that you used.

If the next page of that paper describes how to apply an FEM to a BVP where one "boundary" is infinity, then I need to see that. If it simply involves substituting some finite upper limit of eta for infinity, then this system can be handled by Maple's built-in BVP solver. Indeed, there have been a great many Questions here on MaplePrimes over the past few years about this and closely related ODE BVPs.

@jefryyhalim 

You can get more waves by extending L. I have no idea whether this is physically valid or makes sense in your model; I'm merely describing what's mathematically and computationally possible.

We extend the BVP solution to use x > L. This is done by

1. solving the BVP;
2. evaluating that solution at x=0 to get ICs (initial conditions), and parameter values (the eigenvalue sigma);
3. using those ICs to resolve the system as an IVP (initial value problem);
4. IVP solutions are not bounded to x= 0...L.

Here is the code needed to do this. This starts at the point where you've defined all the ODEs and BCs:

DepVars:= [Vup1,Vp1,Vp2,Vlp1,Vlp2](x):

#Solve BVP (without output= listprocedure):
dsolBVP:= dsolve({deq||(2..6)} union {eq||(1..21)}, numeric, maxmesh= 2^14):

#Evaluate all functions and relevant derivatives at 0. 
#Use D-form to restate them as initial conditions.
#The remove(evalb, ...) is to remove the "x=0" because it's not an IC.
ICs:= remove(evalb, subs(x=0,  convert(dsolBVP(0), D))):

#Resolve same ODEs as an IVP (without output= listprocedure).
#The eval(..., ICs) is to set parameter sigma. 
#The select(hastype, ...) is to remove "sigma= ..." from ICs.
dsolIVP:= dsolve(
   eval({deq||(2..6)}, ICs) union {select(hastype, ICs, function)[]}, 
   numeric
):

#Use side-by-side plots to show equal functions:
plots:-display(
   `<|>`(seq(
      plots:-odeplot(
         dsolIVP, [x,f], x= 0..3*L, view= [0..3*L, -3.2..3.2],
         labels= [x,f], labeldirections= [HORIZONTAL,VERTICAL]
      ), 
      f= DepVars
    ))
);

plots:-odeplot(dsolIVP, `[]`~(x, DepVars), x= 0..3*L, legend= DepVars, size= [1500,500]);

The plots reveal some redundancy that is striking to me: Vlp1(x) and Vlp2(x) are identical, and Vp1(x) and -Vp2(x) are identical. My guess is that this indicates that there's something important missing in your model; however, that's just a guess as I have no knowledge of the physics.

The full worksheet with displayed plots is here:  Longitudinal_ODE_IVP.mw

@jefryyhalim Your sigma is not 0; it merely appears to be 0 at the scale of the plot. To get a more-prescise numeric value, do

A(0);
      0.661316346230992e-2

Does that value make sense in the context of your problem?

@jefryyhalim Regarding the removal of small imaginary parts, I think that at least one of us has had a fundamental misunderstanding of what the other meant. I meant that they could be removed after they had been generated.  So my technique can't speed up the process by which they are generated. That's the part that's taking 5-10 minutes, right?

Regarding the numeric BVP solver, Rouben has provided an example above in the time since you asked. See if you can get that working. The BVP solver in this case is nearly instantaneous, and it doesn't give solutions with any non-real numbers. I think that it'll ultimately provide a more satisfactory overall solution to your problem than will symbolic solution. If you get any error messages from it, do not despair. There are many options for error control that we'll be happy to help you with. I used the option maxmesh= 2^14 for your problem (which is probably a lot more than is needed, but the solution is still nearly instantaneous, so I see no reason to not do it).

@Earl You can also abbreviate Transpose(U) as U^%T, without needing to load LinearAlgebra.

@Rouben Rostamian  This is a Reply to your previous two Replies. In the first, you said "it returns nothing, indicating that there is something wrong with your system." Symbolically solving for the integration constants in this system involves solving (algebraically) a massive system of transcendental equations that is generated by the BCs at the "other" endpoint, the one where the independent variable is not 0. The chances that such a system can be algebraically solved are virtually nil. I wouldn't call this "something wrong with your system." This same thing happens when using symbolic dsolve on all but the most-trivial BVP systems, so it's generally not even worth trying for a symbolic particular solution. You can fsolve for the coefficients after getting the symbolic general solution. Sometimes giving dsolve the option convert_to_exact= false combined with using decimals in the input system helps. 

In the second Reply, you suggest using dsolve(..., numeric) and specifying BCs at three values of the independent variable: 0, L/2, and L. The numeric BVP solver will immediately reject this: It strictly only allows BCs at the endpoints. Compare:

sys:= {diff(y(x), x$3), y(-1) = 1, y(0) = 0, y(1) = (1)}:
dsolve(sys);
dsolve(sys, numeric);

@Christian Wolinski Yes, I thought that you were thinking along the lines of atomic. That's slightly different. Rather than invariance under op, it's invariance under mapping of the identity function, so it's equivalent to

`type/Atomic`:= x-> evalb(x = map(x-> x, x));

I think that everything that's invariant under op is a "fundamental" type, by which I mean a structure that has a direct internal representation (and hence a dagtag  (see dagtag under ?kernelopts)). Of those that are first-class objects (look that up in Wikipedia), by which I means assignable to a variable at the Maple-language level, the only ones that I can think of that are invariant under op are integer, string, and symbol. So, you could use

`type/InvariantUnderOp`:= {integer, string, symbol};

However, this is less robust than the procedure-based methods that I advocated before because if a new fundamental type is added to Maple, then the categorical typespec immediately above will need to be checked and perhaps rewritten.

@Christian Wolinski Your objection to satisfies seems incredibly naive considering the relatively high level of skill in symbolic computation that you've otherwise shown here on MaplePrimes. What do you think "Maple types" are? Were they carved into stone tablets on Mt Sinai and given to Moses? No, they are all procedures that evaluate boolean conditions. Some of those procedures are builtin; most are explicit Maple library procedures named `type/...`. So, if you want it to be more like an official Maple type, just do

`type/InvariantUnderOp`:= x-> evalb(x=op(x));

and install it in the library. 

@jefryyhalim 

@@ is the functional iteration operator. Defining recursively,

(f@g)(x) means f(g(x)) (i.e., functional composition),
(f@@n)(x) (for n::integer, n > 1) means (f@(f@@(n-1))(x),
(f@@1)(x) means f(x), 
(f@@0)(x) means x (i.e., the identity function).

In some limited contexts (i.e., when Maple has been taught what the inverse of f is), 
(f@@n)(x) (for n::integer, n < 0) simplifies to (g@@(-n))(x) where g is the known inverse function of f.

The output of a function may itself be a function. These are sometimes called functional operators in the mathematical literature, but they're so easy to construct in Maple that a special name is not usually used for them. If F is a functional operator, then F(x)(y) means (F(x))(y), (i.e., the application of a function to its argument(s) is itself a left-associative operator). Such is the case with D, the differential operator. So,

D(f)(x) (for x::name) is essentially the same as diff(f(x), x), and
D(f)(a) (for a:algebraic, a::Not(name)) is essentially the same as eval(D(f)('x'), 'x'= a) (or f'(a)). So, for example, (D@@2)(f)(0) means (D(D(f)))(0) means f''(0).

Partial derivatives can be expressed in the functional D form also by applying an index to the D. For example, 
D[2](f)(x,y) means the first derivative of f with respect to its 2^nd-position argument (y in this case), which is essentially the same as diff(f(x,y), y), if x and y are distinct names. This indexing of D can be iterated: Defining recursively,
D[m,n](f) means (D[m](D[n]))(f), which is a second derivative (even if m=n). This indexed form of D can also be used for a univariate function, so D[1,1](f)(0) is another way to express f''(0).
 

@acer I'd expect it to resolve under the integration because the limits of integration contain the necessary information. Indeed, regardless of the presence or absence of assumptions, simplify applied to an unevaluated integral will use assumptions based on the limits. So:

restart:
Psi00:= exp(-r^2/2)/sqrt(Pi):
int(int(conjugate(Psi00*r^(1/4))*Psi00*r^(1/4)*r, r= 0..infinity), phi= 0..2*Pi);
#Above returns only partially evaluated.
simplify(%);