@bliengme 

All Matrix/Vector arithmetic can be done with the familiar binary infix symbol operators: +, -, *, ., /, ^. (Note the distinction between commutative (or scalar) multiplication, *, and noncommutative Matrix/Vector multiplication, `.`.) There is no need to use any package (such as LinearAlgebra) to define Matrices and Vectors or to do arithmetic on them.

Your Reply makes me think that you're now wondering So what's the purpose of multiply and Multiply? 

The command multiply---as well as the entire linalg package that it comes from---are deprecated. That means that they are no longer intended to be used and they have been replaced with better commands. The reason that they still exist in Maple is so that old code that uses them will continue to run.

The LinearAlgebra package has numerous commands with the spelled-out names of arithmetic operations like Multiply. These commands are indeed for performing the arithmetic operations indicated by their names. The only reason to use such a command is if you want to use options, which cannot be included when using a binary infix symbol (because there's literally no place to type the option). The options are described on the commands' help pages. For example, the most-commonly-used option is (I expect) inplace, which means that the result overwrites one of the operands. There is no need to use options unless you're trying to write highly optimized code. Thus there's no reason for a beginning Maple user to pay any attention to the commands with names like Multiply.

Suggestion to Maplesoft: There are too many top-level commands in LinearAlgebra, which makes it unwieldy for even an experienced user to find what they're looking for. All of these not-commonly-used Matrix/Vector arithmetic commands should be bundled into a subpackage.

@Carl Love  So, now that I've read more of the package's help pages, I submit the following replacement for my procedure E__n in the Answer above. I think that it combines all of the good points of Tom's, Acer's, and my own approaches discussed in this thread so far:

restart:
ScientificConstants:-AddConstant(
   `Marc's constant`, symbol= 'E__n', 'derive'= 'h' * 'c' * 'R'[infinity]
);
E__n:= n-> evalf(ScientificConstants:-Constant('E__n', 'units')/(-n^2));

(The constant just needs a better name than `Marc's constant`, which Marc can decide on.)

Note that this procedure tacitly uses the constants to the full precision that they are known modulo the value of Digits at execution time.

@tomleslie You're right, Tom: Introducing evalindets was a bit of overkill, and I was hesitant to do it. Your code still accomplishes my top goals:

• elimination of magic numbers,
• replacing them with widely recognized symbols,
• using them to (at least) the precision that they are known,
• simplification of the units.

If I had had more familiarity with the package, I would've used Constant(..., 'units'), as you did. So I read more of the package's help and wrote a much simpler procedure (posted as a Reply to my Answer) that I think accomplishes all of our goals discussed in this thread.

@acer Thanks, Acer. I considered writing it such that the execution-time value of Digits rather than the definition-time value would be used, but I decided in favor of the (slightly) better execution speed of using the definition-time value. But I certainly see the benefits of doing it your way.

This was my first usage of ScientificConstants, and I found the format that the information is returned is awkward. The return value of GetConstant should be type set(name=anything) or something similar, and the units should be real Units, not simply products of symbols.

It looks like the package has a facility for adding constants, and the constant h*c*R could be added as a "derived constant" (using the terminology from the package's help page).

I figured out the meaning of the constants in the OP's original E__n by Googling them in the form "1.097 x 10^7", etc. I guessed that the 2nd procedure, eV, was simply a units conversion factor, so I simply checked that in Maple.

@waseem Thank you for the pages. So, they just use infinity = 8 for the upper boundary point of eta. In that case, the system can be easily solved with Maple's built-in BVP solver. As I said, this and related BVP systems dealing with nanofluids have been discussed and solved numerous times here on MaplePrimes.

Can you post the whole paper? Or just give me reasonable default values for all parameters?

Before you go through a lot of effort to produce this Fortran code, you should exhaust the possibilities of solving the ODEs in Maple first. Then, if that's not fast enough, do the Fortran. Since you've never posted to MaplePrimes before, it seems unlikely that you've exhausted the possibilities for solution entirely within Maple. We have several regular contributors who are expert in efficient numerical solution of ODE systems. 

Consider these two flowcharts for your project:

Method A:

1. Write some crappy Maple code.
2. Translate it to "optimized" Fortran.
3. Run Fortran code.
4. Does it need to be modified? If yes, go back to 1.

Method B:

1. Write some good Maple code.
2. Run it.
3. Does it need to be modified? If yes, go back to 1.
4. Optimize Maple code. Go to 2.
5. Is it fast enough? If yes, you're done; if no, go to 6.
6. Translate to optimized Fortran.

Method B requires much less effort on the part of the programmer and will produce faster code.

@Kitonum Your plots B and C are transposed. In other words, the xs and ys are exchanged with each other. It is only because of the coincidence that x and y are relatively close that this error was not more obvious.

That x = y - sqrt(y-1) and x = y + sqrt(y-1) are better approximations as y -> infinity can be seen by taking the result of your solve command immediately above and dividing the numerator and denominator through by y^2. Of course, this means that the terms of the radicand are divided by y^4. Then the asymptotic function is obtained by discarding any remaining terms that have any power of y in the denominator.

This does not mean that the result that you got from asympt is wrong, because asymptotic terms are not unique. For example, clearly sqrt(y) and sqrt(y-1) are asymptotic to each other because the infinite limit of their ratio is 1.

@Kitonum Closer asymptotic curves are x = y+sqrt(y-1) and x = y-sqrt(y-1).

@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