Here's how to do it.
 

Download 3dModelFit.mw

The following is a variation of VV's procedure that uses the low-level operations of LinearAlgebra:-Modular for efficiency. I believe that it'll run about 2-1/2 times faster---not a huge improvement, but perhaps worthwhile. Please test it on a matrix known to be MDS.

IsMDS:= proc(
    A::Matrix(square),
    r::And(posint, satisfies(r-> irem(upperbound(A,1), r) = 0))
)
uses LAM= LinearAlgebra:-Modular;
local
    n:= iquo(upperbound(A,1), r), 
    rng:= t-> (t-1)*r+1..t*r,
    M:= LAM:-Mod(2, A, float[8]),
    k, kr, P, ip, jp, det
;
    for k to n do
        kr:= k*r;
        P:= combinat:-choose(n,k);
        for ip in P do 
            for jp in P do
                LAM:-RowReduce(
                    2, M[rng~(ip), rng~(jp)], kr$3, 'det', 0$4, false
                );
                if det=0 then return false fi
            od
        od
    od;
    true
end proc
:    

My try:

restart:
d:= proc(m::integer)
option remember;
local i;
    -expand(add(xi[i]*thisproc(m-i+1), i= 2..m))
end proc:
d(1):= 1:

DD:= m-> rtable((1..m)$2, (i,j)-> d(j-i+1), subtype= Matrix): 
DD(4);

The command gcd is meant for polynomials. As polynomials, a and b have no common factors, so their GCD is 1. The command for the GCD of integers is igcd.

You make a, b, and x the parameters of a procedure and pass to them two expressions and a variable name. To do RK4, we also need an initial condition, a final value of x, and a number of steps. Like this:

MyRK4:= proc(
    a::algebraic, b::algebraic, x::name,
    x0::realcons, y0::realcons, xf::realcons, n::posint
)
local y, h:= (xf-x0)/n, k;
    dsolve(
        {diff(y(x),x) = 10*(y(x)-a)/b, y(x0) = y0}, numeric,
        'method'= 'classical'['rk4'], 'stepsize'= h,
        'output'= Array(1..n, k-> x0+k*h)
    )[2][1]
end proc
:
MyRK4(ln(x), x, x, 1, 1, 2, 10);

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.