@zphaze I think that you're failing to see the significance of acer's 2nd example. He did a restart, so this example is not using any package at all. That's not quite what I'd call "switching packages".

@adel-00 I've made numerous Posts and Answers here over the years showing 3d plots of parameterized dsolve solutions where one or both of the.plot's independent variables are parameters of the ODE system. I prefer to use plot3d for this (as shown by mmcdara (vote up)) rather than surfdata because plot3d lets you view the parameter as a continuous variable and because using surfdata feels too much like mucking around in the low-level details of plot structure. 

@nm It is well known that some Maple results are session dependent. This particularly affects the order of results, but it may also affect the form. This can mean that solution pathways through algorithms may depend on random factors or factors out of the user's control. Are you sure that that is not what you're experiencing, with the warnlevel thing just being an unrelated coincidental factor?

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