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

@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.