@anton_dys I don't know whether the following will more help you in the short term or confuse things in when going awry for more complicated examples, but...

Your given 3D plot generated using the plot3d command can be converted into a 2D density plot (like the surfdata or densityplot commands generate) without having to do recomputation. And the resulting 2D PLOT structure can even be further adjusted using plots:-display with its size option. And since the original call is to plot3d rather than surfdata there is no need to discretize the expression like sin(x*y) to a Matrix or 2d Array data structure as I did above.

Not all 3D plots can be converted in this way. It works for yours primarily because the data is in a GRID substructure and the surface was shaded by a custon shading scheme.

Having said all that, I still think it might be more prudent for any user who is not conversant with the subtleties of Maple's plot structures to instead rely on the surfdata approach I gave in my original answer (until such time, say, as the densityplot command might respect the colorscheme option differently).

restart;

P3d:=plot3d(sin(x*y), x=-5..5, y=-5..5, view=-10..10,
            colorscheme=["zgradient",["blue","green","yellow","red"],
                         zrange=-5..5],
            style=surface, lightmodel=none, orientation=[90,0,0], axes=box):

P2d:=':-PLOT'(remove(type,[op(P3d)],
              ':-specfunc(anything,{:-AXESLABELS,:-LIGHTMODEL,:-ORIENTATION})')[]):

P2d;

plots:-display(P2d, size=[600, 300]);

@anton_dys I have a pending post on 3D plots and size and aspect ratios. 

For 3D plots the size of the (projected, flat) display window, and the aspect ratios of the 3 axes should be controllable by different means. It would be a terrible mistake if some gui developer conflated them. This special case at hand, with an orientation that happens to line up with the axes, must NOT be used to characterize the wider need.

@Joe Riel A note of caution about pre tags: when pasting plaintext code within pre tags I have to replace lessthan symbols with ampersand-lt-semicolon entities, so that the editor doesn't get confused.

I toggle to html mode in the editor (using its icon) and then paste into a pre tag.

acer

@Age I had indeed noticed the problem above rho=50 but not looked more closely because the solution for your IVP at hand hand rho(r) stay below rho=8.

I can take another look, but it may be a few days before I get a spare evening.

I ought to have asked earlier: what kind of accuracy would be acceptable in the solution to the IVP? Are you hoping for a solution with merely the same general qualitative feel, or if something numerically accurate then how accurate? I ask because you have numeric integration and rootfinding within differentiation (possibly with a numapprox layer) within numeric ODE solving. Often the working precision and accuracy requirements will need to get progressively tighter as one goes deeper in the level of approximation, to ensure that the highest result (the IVP solution) is not nonsense. Adjust wrongly at any step and final or iterim results can become rubbish or take forever to compute. I have been trying to balance all that so far, and also to attain reasonable speed performance, but it's hard not to second-guess everything unless I have an idea of the overall target.

There are some things that I didn't understand early in the sheet (the use of rho=250 and the 125^2 near where `func` was defined, etc).

I may rewrite the whole thing in 1D notation in a Worksheet because the irritations of 2D Input are just too much when I am trying to think hard.

@Age Various edits and a few alternative approaches.

The plots look less ragged in Maple itself. I ran this in 64bit Maple 2015.2 on an i5-6600K.

For some reason the mapleprime's inlining mechanism doesn't seem to be showing the timing printouts from the wrapping CodeTools:-Usage calls (but you can see them in the document).

Download ds_help_1.mw

I expect that we could sort out that derivative for you, too, so that you don't have to rely on such a crude approximation.

And we may be able to speed up or improve the integration too -- my dsolve hasn't completed even yet, and it's not clear that the numeric integration is succeeding for all rho(r). Give us a little time.

Perhaps next question you could go straight to the real document & issues, thanks.

@Age Uncomment the line in ODEproc which allows it to return unevaluated for symbolic argument. But correct it with respect to Greek letter rho vs English letter p.

Also correct your initial conditions (ICs) with respect to those names. Your sheet has English p for the ICs but uses typeset Greek rho in the DE set `Eq`.

Maybe you should consider using a Worksheet with 1D Maple Notation for input rather than a Document with 2D Math for input, if this kind of thing gets in your way.

In which platform and Maple version (major.minor number) are you seeing this unwanted switch from Text mode to 2D Math mode?

acer

I could add that I was able to get an exact result for the integral in question, and but even after quite a bit of successful simplification it contained large exact coefficients which make it prone to round-off error unless Digits is set too high to allow the bvp solver to compute quickly. Also, float evaluation of the exact integration result (simplified ot not) was giving a 0.0*I artefact that was awkward.

I also looked briefly at a series solution for dealing with the integral term (trying for both the integrand as well as simplified exact integration result mentioned). I don't like this approach (in general and here in particular) because of concerns with accuracy of its approximation, especially without knowing the plausible range of f3(x) of a dsolve solution.

@vv The first of these I see wrong back to Maple 9.

is( y<=0 ) assuming Or(y<0, y=0, y>0);
                                                  true

The second I see as OK in Maple 14.01 but wrong in Maple 15.00.

is( y<=0 ) assuming Or(Or(y<0, y=0), y>0);           
                                                  false

kernelopts(version);
                        Maple 14.01, X86 64 LINUX, Oct 15 2010, Build ID 535952

and,

is( y<=0 ) assuming Or(Or(y<0, y=0), y>0);
                                                  true

kernelopts(version);
                        Maple 15.00, X86 64 LINUX, Mar 20 2011, Build ID 582305

I'll submit a bug report for both.

@AndrewG One difference is that you can assign the result from RandomVariable to a name, and then use that name in multiple instances (including one or more calls to Sample). But sometimes you might not wish or need to make that much effort. You might want to create a sample as a one-off operation, in which case a call to Sample(Normal(...)) is a more convenient syntax.

@vv My results could be "fortuitously" predicated upon this behavior,

restart:

is( y<=0 ) assuming Or(y<0,y>0,y=0);

                                                           false

is( y<=0 ) assuming Or(y<0,y=0,y>0);

                                                            true                                     

simplify((y^2)^(1/2)) assuming Or(y<0,y>0,y=0);
                                                           | y |

simplify((y^2)^(1/2)) assuming Or(y<0,y=0,y>0); # !!

                                                             -y

@vv The latter cases, where it gives up at csgn, or refuses to combine a product or two radicals, may be indicate one or more parts of the problem. (It's just a guess, but the former case might relate to something baulking on _Envsignum0...)

Download csgnhmm.mw

Did you try it using Int instead of int? That should avoid the symbolic integration and allow fsolve to use evalf(Int(...)) which is numeric integration.

You might need to use the digits and epsilon options in the call to Int, to ensure that the numeric integral is finely enough resolved for fsolve to satisfy itself about accuracy of the root. Quite often one needs the result from evalf(Int(...)) to be at least one or more decimal digits more accurate than Digits, for this use of fsolve.

acer