@MapleMathMatt Your first method, using indices, is very inefficient compared to hasindex or assigned, the lookup time through a list being proportional to the list length (O(n)), and also the creation of the list is O(n). One of the most important things to know for efficient Maple programming is that table index lookup uses nearly constant time (O(1)) regardless of table size (and that constant is very small). It's unfortunate that documentation of this is buried in a terse and densely packed appendix of the Maple Programming Guide.

If you had a situation where index memberships checks occurred significantly more frequently than index additions or deletions, then it might be worthwhile to use the set of indices for the lookups. The ratio of checks to changes would need to be at least o(n) (little-o not big-O) for this to be worthwhile. 

What is `Y1"`? Is it a function that you've already defined?

My guess is that the difference is not in dsolve itself but in simplify (or something related to that). 

This Reply is only to point out an issue with your terminology. I don't at the moment have any solution advice (mostly because I'm on my phone, not Maple). 

Your problem is a boundary-value problem (BVP), not an initial-value problem (IVP). While that distinction may seem like a triviality, BVPs are often much more difficult to solve. 

Your last line of code contains a variable j which isn't used elsewhere. What is it? Is it the same as the imaginary unit I?

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