I haven't (yet) tried to assess the correctness of your argument. This is perhaps a known result. Specifically, it's known that if |p - q| < n^(1/3), then n is factorable in polynomial time.

@Preben Alsholm Thank you. I've made the correction to the Answer.

@J F Ogilvie  Yes, I agree about the old Spread package. I'd used it several times, and it seemed much easier to use programmatically than anything somewhat equivalent currently available through DocumentTools. My guess as to why it became deprecated is that it didn't integrate well with the extensive overall development of interactive embedded worksheet components, it being one of the earliest; I think that only plots are older, and their interactivity is minimal. 

It doesn't matter to me whether the GUI aspects of the old Spread are maintained. I'd be happy having any two-way correspondence between Matrices and an interactive GUI component that supports scrollable arbitrary width.

@Lyssaloo There are at least 3 different direction fields that could be plotted in your case: Er vs t, On vs t, and Er vs Or. DEplot just needs to know which to plot.

@mmcdara I wasn't referring to Maple's piecewise command. I mean that all functions of a real variable that are used in Maple (or any CAS) in any practical way are continuous or piecewise continuous in the mathematical sense. Thus the set of discontinuities is countable and nowhere dense, and it's theoretically possible that discont could report them all so that we can then assume that the function is continuous on the remaining intervals.

@adel-00 

r1, which is the end result of all our work in this thread, is nonetheless just a simple multivariate polynomial with float coefficients. Thus any further process that starts with r1 ​​​and ends with results that you question is just about simple multivariate polynomials and has nothing to do with the foregoing process of double integration, etc.

I just want you to be clear that this RootOf thing is not caused by any of the advice that you've received so far in this thread. 

Generally, polynomial solutions returned as RootOf expressions with float coefficients can be further simplified with allvalues. I don't know about this case; I haven't looked at it in Maple. If there are several unspecified variables multiplied together, it may be impossible to proceed without giving them numeric values. There's also a chance that giving solve the options parametric and real will provide a useful answer. But it may also give an overwhelmingly lengthy result. 

@adel-00 In some cases, this being one of them, the expression needs to be preprocessed by collect(..., distributed) before being passed to coeffs. Since it's not clearly documented when this extra step is needed, and since it can be more computation intensive than the coeffs itself, I'll use the following method to do the collect only if the coeffs at first fails. Change procedure Coeffs to

Coeffs:= proc(e::algebraic, V::{list,set}(name))
local T, C;
    try C:= coeffs(e, V, T)
    catch "invalid argument":
        C:= coeffs(collect(e, V, ':-distributed'), V, T)
    end try;
    table([T]=~[C])
end proc
:

Then following through with the rest of the code, we get

r1 := 0.6218523429*a*b*alpha1 - 3.902805269*a*alpha2*b^3 - 3.825073736*b^2*beta1 + 
    26.09719472*a^4*beta2 - 11.78614736*b*a + 26.09719472*b^4*beta2 + 
    48.21385264*a^2*b^2*beta2 - 3.825073736*a^2*beta1 - 3.902805269*a^3*b*alpha2

Of course, this process can only factor the symbolic constants out of the integrals if it's mathematically possible to do that.

@dim____ 

Then the hint that you were given is a bit off the mark. While mul can be made to work for this in the way shown by nm, the commands product or Product (its inert form) are far more appropriate for a situation where the number of factors being multiplied is not exactly specified.

This formal analogy of Maple commands is exact:

(product : mul) :: ​​​​​(sum : add)

The help pages for mul and add are the same page.

@adel-00 

Regarding double and multiple integrals: Yes, as a syntactic convenience, double (and higher order) iterated integrals can be specified with a single invocation of int with the second argument being a list(name= range), as shown. This is completely analogous to the following two equivalent ways of specifying a mixed 2nd partial derivative:

diff(f(x,y), [x, y]);
diff(diff(f(x,y), x), y);

Regarding additional symbolic coefficients: Yes, the new coefficients can be accomodated with a very simple change to my code: Change the 2nd argument of add from v= V to v= [indices](C, nolist). (Of course, the new coefficients must be included in list V also.) Indeed, this change will work better for your earlier case anyway, so I'm changing the Answer to this new form. Note that my original Answer had a mistake in that it didn't include the term that has no symbolic coefficient.

@Kitonum 

My command max(abs~(F-G)) is not equivalent to the matrix infinity-norm of F-G. The infinity-norm of a matrix is the maximum of the vector 1-norms of its rows. My command gives the infinity-norm of F-G viewed as a vector.

While I do realize that a variety of norms could be used to measure the "distance" between matrices, I don't see how any matrix norm could represent "the largest distance between the two curves [surfaces]." In particular, note the word largest, which also appears in the thread title. Likewise, I don't see how any vector norm other than the infinity-norm would be appropriate for measuring the largest distance.

Observe the difference in the following two infinity-norm computations. In the second, Alias is used to recast the matrix as a vector.

A:= LinearAlgebra:-RandomMatrix(2);
                          [-63  30]
                     A := [       ]
                          [-26  10]

LinearAlgebra:-Norm(A, infinity);
                               93

LinearAlgebra:-Norm(ArrayTools:-Alias(A, [numelems(A)]), infinity);
                               63

@mmcdara This type of error occurs in 2D Input, regardless of whether it's document or worksheet mode. (I think that you may think that 2D Input and document mode are the same thing.)

Anyway, vote up.

To the OP: When using 2D Input, if you put a space between a function's name and the left parenthesis, then the function name is treated as an ordinary (scalar) variable. Conversely, if you have a variable coefficient immediately followed by a left parenthesis, then it'll be considered a function rather than a coefficient. I think that these two things are the most common errors made in 2D Input.

The indets command can be a great help in finding these errors. 

@Jaqr The list L has 3 elements, but you're indexing it as if it had 7 elements.

@Assinat Using 1D input, ^+ means transpose. In 2D Input, use ^%T instead. The latter form will work regardless of input mode.

AFAIK, assuming is something that only makes sense for variables or expressions containing variables; you're trying it on a constant. IMO ideally assuming a::posint should return an error if a is constant. I don't know why it doesn't. @ecterrab Edgardo, do you have a comment or opinion?

@Kitonum @vv

Yes, my solution is only valid when there are two real solutions to h=f in a..b. This was my oversight. The critical p intervals can be determined using the parametric option to solve, among other ways.

@Katatonia : Yes the function g can be expressed with piecewise. ​​​​​This makes the raw result of the integration far more complicated, but it simplifies to the same thing.