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

@666 jvbasha If I execute your worsheet top-down, I get "Error, (in evalf/int) invalid arguments". This suggests that the output that you show is not the direct result of the code that you show. That could possibly be the result of executing things out of order.

@mmcdara I included the column for your "reduced uncertainty", but I just called it "z-score" because I think that's what it's more widely called (at least in English). I also kept the p-value. I changed the **** flag to a variable-length number of stars (0 - 5) depending on how low p is. I added an option to test either the left or right side (with the default being the right), which is important for the nonsymmetric distributions. You can see all this in my new Answer.

@epostma Hi Erik,

Please see my Answer below where I've tested every stock continuous distribution for this sampling irregularity. The ones that have it are ChiSquare, Error, FRatio, InverseGaussian, LogNormal, and Normal. The Exponential distribution does not show this anomaly, so I guess that the problem is in the "fallback algorithm for the tail", which isn't needed for the Exponential.

Let me know if you have any questions about my testing method.