Use Expand instead of Normal.

Why Expand works and Normal doesn't: Suppose you remove the Normal and just do the Evals. Then you have substituted values for all the xs and ys. The only variable left is q. But q isn't really a variable; it's a constant field element. So the expression isn't really a rational function (except in a degenerate sense); it's just one field element divided by another and is thus just a field element itself. The division of (higher-degree) field elements can be done by Expand, but apparently not by Normal. Personally, I think that this is cruddy design, and Normal should work for the degenerate cases also, but that's the way it is.

Note to the other correspondents: This assumes that p has been assigned an explicit prime value, and q has been aliased to the RootOf of a polynomial irreducible mod p. Having Answered many of the OP's Questions over the past few weeks, I am certain that they had done these two preliminary steps before executing the code shown in the Question.

Note to the OP: People will be more likely to understand your Questions if you post all the code necessary to reproduce your results or error messages.

The operator precedence of mod is higher than that of =. Thus Q2 = P2 mod p is equivalent to Q2 = (P2 mod p). To get it to do what you want, use (Q2 = P2) mod p (as suggested by Christian). The complete list of operator precedence is on help page ?operators,precedence.

I need to see EAdd to be sure, but my guess is that there is some exceptional case where it does not return a pair. Could it be the point at infinity? 

To have a procedure remember all its past arguments and return values, include

option remember; 

in the declarations section. To access this remembered information, use

op(4, eval(P))

where P is the procedure's name. The information is not stored in the order that it was generated.

@nb99 You wrote:

• I'm somewhat confused as to the machinery of asympt.

You can look at the code via showstat(asympt). It's fairly simple for a Maple library procedure of "higher" math. Basically, it computes the series about z=0 of f(1/z), but I do think that it also incorporates some special knowledge for some special functions. I don't know where that knowledge is stored, but no doubt that knowledge-base can be found eventually by digging through the code.

• Indeed, for your case it works and generates a converging power series. However, if one replaces in the argument of HeunC 1/z -> z, for example, it fails.

Yes, I consider that I "cheated" by using 1/z, especially considering the previous paragraph above. I just wanted to show that some asymptotic series for HeunC was possible while placing very minimal restrictions on the parameters.

• The naive explanation would be that there isn't a converging power series for that case.

A superficial glance at the series leads me to think that that's the case: Replacing the parameters with with numbers, the coefficients do not appear to be going down in magnitude at a sufficient rate to have absolute convergence for |z| > 1.

• However, looking at the documentation of asympt, it seems to be able to generate more complicated asymptotic forms (e.g. a cosine/sine with decreasing envelope for a Bessel function). Can't it find, say, an exponential dependence?

Yes. Probably the most well-known example would be asympt(n!, n). But note how the result is presented with n^n*exp(-n) factored out of a SERIES object (note where the O term is). This is a small clue to how the machinery behind asympt works.

• A second hypothesis is that it might be problematic to use asympt when the argument of HeunC is greater than 1. Could that be the case?

That may be the case as far as Maple's asympt is concerned. But you say that certain asymptotic expansions are known? And for an arbitrary number of terms in the series? If so, I'd like to get Maple to spit out at least one of those known cases. Note that all of Maple's power-series-related commands^[*1] return an arbitrary-yet-still-finite number of terms. So, please show me one simple such case from "the literature".

There is a package for power series that's in beta testing (or perhaps it's gamma testing?), but is packaged with your stock Maple nonetheless. It's called MultiSeries. It has its own asympt command. (See ?MultiSeries.) That package has hooks for adding knowledge of new functions. Currently, it has no knowledge of HeunC (see MultiSeries:-FunctionSupported(HeunC)). Perhaps we could add sufficient knowledge under some severe restrictions on the parameters and get a useful series.

[*1] The one exception that I'm aware of is convert(..., FormalPowerSeries), which returns true infinite series. But the cases where that works are very limited.

An alternative to stopat is stopwhen, and I find it much easier to use. Using it, you can explicitly set breakpoints in your code, and they won't be keyed to those ephemeral "statement numbers" (that may change every time you alter the procedure). Here's how. Pick a global name. I chose DBG_BP (for debug breakpoint). In your initialization file, put the commands:

DBG_BP:= 0;  stopwhen(DBG_BP);

Now, any place in your code that you want to set a breakpoint, include the command

++:-DBG_BP;  #That's Maple 2019 syntax

Whenever that command is executed, the debugger will pop up.

You can get even more flexibility by putting the breakpoint in a userinfo statement, like this:

userinfo(1, 'DEBUG', ++:-DBG_BP);

The key number 1 can be any positive integer of your choosing, and the keyword DEBUG can be any of your choosing. When you want to turn on debugging, give this command before executing the code:

:-infolevel['DEBUG']:= 1: #or any number >= the key number

To turn debugging off, 

:-infolevel['DEBUG']:= 0: #or any number < the key number

You can have multiple key numbers for different "intensity levels" of debugging.

The values of any variables can displayed with userinfo statements, even plots.

I can't tell which, if any, parts of your computation are meant to be inert (as in Diff rather than diff and Sum rather than sum). This may be simply due to my lack of sufficient knowledge of document mode. (In worksheet mode, inert operators are gray instead of blue.) I don't think that there's any point in trying to compute the derivative of a sum with respect to one of its limits of summation unless Maple can find a closed form for the sum. So, my suggestion is that you work on the sum separately.

It would be nice of you to post the original system. 

From a quick glance at your results, it appears that there's a repeated eigenvalue, -1, and that would make it impossible to put the result into the form that you stated in your second-to-last line.

You have a system of 5 ODEs with only 4 function variables: f(eta), G(eta), theta(eta), and H(eta). You must remove 1 ODE. Since the first 3 ODEs only contain f(eta) and G(eta), you should remove one of the first 3. Indeed, by differentiating Eq3, it's clear that Eq2 and Eq3 are contradictory. So, one of those 2 must go. 

Ah, I can tell now that you're referring to Maple Companion because of the tessellated background in your Question.

As far as I can tell, Bessel functions are not implemented there. [Edit: They are implemented. See the Replies.] Maybe someone else here knows for sure. However, functions are not limited to what's on the functions keyboard. For example, erf is implemented and can be entered from the alphabetic keyboard. I tested erf(1.).

Regular Maple 2019 can do your integral (using Acer's second interpretation, which seems far more likely than his first).

All functions should have their arguments in parentheses, e.g., ln(x), not ln x.

Try sending your big plots (here big means a lot of data for a plot) to a separate window. Before giving the plot command, use

plotsetup(window);

For example, I just plotted this in a window:

plot3d(sin(x)*cos(y), x= -Pi..Pi, y= -Pi..Pi, grid= [300,300]);

which is more than twice as much data as your matrixplot. The plot appeared instantly, and rotations were fast and totally smooth.

To go back to plotting in the worksheet, use

plotsetup(default);

Also, I've noticed that the speed of all worksheet operations slow as the total size of the worksheet grows. The separate plot windows are not part of the worksheet.

Skip the assume. It's not needed and may interfere with getting a solution (I didn't check). The assumption is effectively stated in the fsolve command:

fsolve(eq, x= 49.32883964..infinity);
                          198.3826229

As a general rule, making assumptions on variables that are being solved for is not a good idea. Making assumptions on other variables that may appear in equations being solved is another matter.
 

There are many ways to do the problem, and CurveFitting:-ArrayInterpolation wouldn't be my first choice. But, since you've already started with it, I'll finish with it:

restart:
V_Ed13:= <257.184, 230.4, 184.3, 138.2, 92.2, 46.1, 0>:
x_12:= <-149, 208, 567, 925, 1283, 1642, 2000>:
V_Rdc:= 153.25:
CurveFitting:-ArrayInterpolation(
   V_Ed13[(P:= sort(V_Ed13, output= permutation))], x_12[P], V_Rdc
);
                        808.125813449024

The sort is needed because ArrayInterpolation wants the independent-variable values in ascending order. The permutation P is to rearrange the x_12s (which are now the dependent variable) exactly the same way as the sort rearranged the V_Ed13s (maintaining the pairwise correspondence). I realize that in this case I only needed to reverse the order of the array with a minus sign, but I wanted to be more general. So, in this case, this gives the same result:

CurveFitting:-ArrayInterpolation(-V_Ed13, x_12, -V_Rdc);

Both of these sort-and-transpose solutions will only work (i.e., produce correct results) if the interpolating function is monotonic. In this case, it clearly is. 

The easiest way is elementwise operators. These are operators whose operands may be containers such as lists but whose action is on the containers' elements rather than on the containers themselves. All elementwise operators end with the tilde character ~.

A:= [1,2,3,4,5,6,7,8,9,10]:
B:= [10,20,30,40,50,60,70,80,90,100]:
E:= 5:

C:= A*~E/~4*~Pi*~B:

See help page ?elementwise. Note that the operands are not required to be containers. The E, 4, and Pi above are not. 

Any user-defined function or existing Maple function can also be made elementwise by appending ~ to its name:

C:= (A,E,B)-> A*E/4*Pi*B:
DD:= C~(A,E,B);

So, in the above the elementwise operator acts much like zip (a command described in the other Answers). But it is more powerful than zip because it works for any number of operands, not just two.

Keyword parameters can be declared and used in a much easier way than as shown in Christian's Answer. You almost had the syntax right, but the curly braces need to go around the entire parameter declaration:

PlotGraph:= proc(
    func::anything,
    {zoom::list({range(realcons), identical(DEFAULT)}):= [DEFAULT$2]}
)
    plot(func, 'view'= zoom, _rest)
end proc
:

Example usage:

PlotGraph(x^2, x= -7..7, zoom= [DEFAULT, 0..7], color= green);

So, note that by using _rest all other passed arguments are also sent to plot regardless of whether they are keyword parameters.

Keyword parameters are used in the body of the procedure exactly like ordinary parameters. The differences are only in how they are declared and passed. Within the procedure, the parameter's value is only the right side of the passed equation; the value doesn't include the parameter's name or the equals sign.

There are several extensive help pages on parameters. You can start reading at ?parameters. That will contain links to the other pages.

Once the 14 equations (5 primaries, 5 solutions, and 4 intermediates) are entered, it's easy to check whether the equations are satisfied. The technique is shown in the worksheet below. Unfortunately, your proposed solution satisfies 4, but not 5, of the equations (assuming I entered them all correctly, which you should check).

It also easy to get a valid solution, but it won't be simplified through intermediate equations such as the one in your original presentation in the PDF.
 

Download eliminate.mw