@HS There are packages for high-speed computations with univariate and bivariate polynomials over finite fields. These are called modp1 and modp2. Most of the functionality of these packages is available through end-user-level commands, for example, mod. But, sometimes you want to work at the nitty-gritty level and need these packages directly. ConvertIn is a command in these packages (and of Maple modules (not algebraic modules) created with GF). I'll agree that the documentation of these is sketchy. See ?modp1 and ?GF. There is a fundamental Maple data structure called zppoly used by these, which is mentioned.in the original Question. 

I don't have an answer for you, but I just wanted to praise you for a very well-asked Question from a first-timer.

I added tag Physics, which is always a good idea for Questions regarding that package. 

@acer Indeed the time difference is disturbing! There's something in both the original code and mine that takes significantly longer in Maple 2016+.

Regarding RAM, I could write it to use less memory at the expense of some time by using an iterator for combinat:-choose. But I'm hoping for an intermediate-sized test case before I do that. 

@Magma There's a bug in older 2D input such that it can't handle some instances of (). You can change that to 0 or anything else. The entries of that table are insignificant; only the indices matter. It's a tactic that I often use: I use () when the entries of a table don't matter. Anything could be used. Most people write NULL.

Do you have an intermediate-sized example, larger than 5x6 but smaller than the one immediately above?

Why do you use overload?

@Magma I'm beginning to work on your procedure now. The first phase is to change the things that are slow due to idiosyncracies of Maple. However, the algorithm itself (as contrasted with simply its Maple implementation) appears to have 5-deep nested loops. I may not be able to overcome that to handle 32x32 input in reasonable time.

@ActiveUser For output, you can use any Unicode symbol. You can even use a Chinese character. Amazingly, this even works in the command-line interface. You just need the Unicode byte values of the character, which you can find on Wikipedia. 

@Kitonum I was supposing 0 points for unanswered questions, so -9 comes from 9 incorrect answers and 1 unanswered. 

@666 jvbasha Sorry, I can't help you. Hopefully someone else can.

@666 jvbasha Sorry, but I don't have sufficient background to begin reading the paper at Eq. 16 and understand it. Perhaps if I had the whole paper, I could.

And why does the plot's caption say "3rd- and 4th-order approximation" whereas its legend says "13th-, 14th-, 15th-order approximation"?

@fatemeh1090 

First, assign the dsolve results to a variable, as

sol:= dsolve(...);

Then, to plot p(t), use

plots:-odeplot(sol, [t, p(t)], t= 0..1);

and do likewise for the other functions. Ordinarily, I'd put the functions together in a single plot, but it's not possible in this case due to extreme differences in the scaling of the vertical axis.

@HS There would need to be some fairly special relationship among the coefficients for the gcd to not be 1. 

@eggs1996 If you have an improved procedure that you'd like further help with, let me know. 

Are you expecting a result other than 1, such as x^3 - e where e is some complicated expression that depends on a, b, and c?

@syhue Yes, my main procedure works for any finite number of dimensions. The extra plotting procedures are for 2-D data.