You might consider a revision to your first (symbolic series dsolve) approach. You could instead search your Y1 list with a binary search. On my machine that reduces its time from 41.1 sec down to just 0.33 sec.

Note also that Y gets formed by merely splitting the range -2..5  by increments of 0.002. So the first element of Y1 which makes exp1>exp2 is a coarse root approximation. A little hammering can make a procedure with which fsolve can find finer resolution roots for your X1 values. On my machine that takes 3.9 sec.

Proc_Error_ac2.mw

If you also adjust your second approach to do a binary search (of Y1) then you can quickly determine that the numeric integration is not succeeding for some avar,bvar values. I'm not sure how that integration issue might be best resolved, or even whether it first occurs before exp1>exp2 is satisfied, etc. But this attachment also shows such problematic values.

Proc_Error_ac.mw

(Could the DE be augmented, instead of doing a separate integration involving the dsolve solution? I didn't really look, to see if that makes sense... Sorry, no time)

A minor adjustment can switch the x- and y- roles in the data, before the interpolation. This flips the densityplots across the x=y line (...one of several ways to do that).

Also, in your fourth construction your final call to plots:-display is wrongly using C11 (from earlier), instead of the new C22 as in that fourth part. I've corrected that. Now the fourth pair agree, with the switch.

Help_4_ac.mw

For example, here is that fourth pair overlaid, with the variable switch before interpolation as well as the C11->C22 correction.

(This code looks like it may have been a modification of this.)

ps. I used Maple 2015.2 for this, since that's your version. In later Maple versions some aspects might be easier or more flexible.

@salim-barzani Another choice for experimenting with plotting options is to use the PlotBuilder (either by right-click on an expression/output or -- my preference -- using the command).

PB_3d_ex.mw

If you put the mouse cursor/focus back on the PlotBuilder's plot then the right-panel gets repopulated with its menus, etc. But you can also switch back-and-forth between that and other work areas in your worksheet.

It's not clear what you mean by "template" above. Are you saying that you want to extract the numeric data generated in the plot structure, and then export that (as a Matrix of numbers, say), to Matlab or other programs? That's trickier for densityplot, but for plot3d you can use the plottools:-getdata and ExportMatrix commands.

For densityplot, one could extract the genetrated numeric Array of shading color values. Let us know if that's what you need, in say your Maple 2021.

You can combine multiple 3D plots by passing them to the plots:-display command.

I'm marking this Question as being for the "Maple 2021" Product, since that's your version (and, for some commands, the version matters). In future, could you please do that yourself? Also, please don't mark your queries about Maple as being for "Mapleprimes" and other products, thanks.

The GUI does not actually use a Maple command, when it adjusts the plot's characteristics via that right-click menu. What the GUI has at hand is the PLOT call structure.

There is some new-ish "re-draw" stuff related to a command that the Maple plotting commands sometimes inject as text into the plotting structure -- eg. PLOT function call, as structure. It can use that now for some zooming. For example you can zoom out of the following plot, and magically the (post Maple 2022) GUI will cause the plot command to be reformed with a widened range, and reexecuted.

P:=plot(sin(1/x), x=0..1, gridlines=true):
indets(P,"input"=anything); # outputs the following:

{"input" = [table([1 = plot, 2 = [sin(1/x)], 3 = (x = 0 .. 1), 4 = (gridlines = true)]), "originalview" = [0.0000192831942484872343 .. 0.997499999997671694, -1. .. 1.], "originalaxesticks" = AXESTICKS(DEFAULT, DEFAULT, _ATTRIBUTE("source" = "mathdefault"))]}

So, that is actually stored in the PLOT function call (data structure) that the Library plot command produces, in modern Maple. And when the rendered, inline plot is zoomed out then the GUI is now able to send a modified reconstruction of the plotting command call, to the kernel. The original structure only had data computed for x=0..1, but after zoom-out the reexecuted call can recalculate for the wider range.

But the GUI's knowledge of how to adjust the call is limited. It knows how to change the range, but not everything. And not all plotting commands store the details in the structure. And, importantly, there's no place for the GUI to hand back the modified command (or even the new output structure) to the user.

A plot in a PlotComponent does allow the user to access the GUI's stored structure. If you right-click the above plot when inside a PlotComponent, and change the color, then you can programmatically extract the modified plot structure from the embedded component (GetProperty). But unfortunately the spiffy zoom-out redraw/reexecute thing doesn't seem to work in a component. And the GUI lacks knowledge about how to itself adjust the relevant stored command data, in most cases.

That leaves the PlotBuilder. You can issue,
   PlotBuilder(lhs(eq));
and, in the right-panel, modify the range and the gridline choices, etc, and toggle on the "Show Command" button. That can be a handy way to explore various option choices, and see/copy the underlying command. Unfortunately it doesn't yet support discrete data (eg. your list of numeric root values, for which pointplot could be a choice).

@salim-barzani In the attached variant, the lastest parameter values used in the exploration are available (outside the exploration). They are assigned in a list to the name last.

graph_in_simple_shape_accc.mw

You could, of course, subsequently programmatically substitute those values into your M expression, if you also need that.

Alternatively, the following variant has the latest instance of M (substituted) assigned to the name last.

graph_in_simple_shape_acccB.mw

Download r01.mw

You can supply 

    output=L

to the,

   CodeTools:-Usage

command call. The list L can contain items like, value, realtime, bytesalloc, etc.

That makes it return an expression sequence of all those output choices. So you can get your hands on the timing/mem results programatically, as additional returned results.

You can then store items in that returned expression sequence, as you see fit. And your computation's own result is included in that sequence, if value is an item in the list L.

That's in contrast to the default behavior of a merely printed message containing the various stats, and only the value as single returned result.

Is this something like what you're after? (I realize that the final construct is not in terms of functional differential operators.)

I've used PDEtools:-undeclare(prime) to get the derivatives displayed in a "subscripted" notation.

Download BE_ac.mw

I realize that this has some flaws; it can't get by discontinuities in the numeric solving, due to the set-up. But perhaps it might inspire you.

(I can't help but feel that I once knew some better way, but right now cannot recall it. I'm sure I'll kick myself, shortly.)

Download Dipole_fields_ac.mw

Is this the kind of thing you're after?

Download test_sol_for_PDE1_ac.mw

You have mixed up the order of the arguments. You've used a different order, in two different places.

You have DOE set up with the columns (apparently) being:
    epsilon, A, Br, N, beta, m, lambda0, lambda1

But you constructed your procedure Fop to treat its arguments as if the order were,
    A, Br, N, beta, epsilon, lambda[0], lambda[1], m

I'm guessing that the wrong one is your construction of Fop, because there you just used a lexicographic order, and ignored how you set up the columns of DOE.

So I suppose that you wanted Fop more like this: Answer_not_same_ac.mw

Such a placeholder is not required. And you don't need to substitute into a diff(...) result.

Instead, you can just use D(W). You can use it as D(W)(u(t)) in your DE for dsolve. And you can use it for D(W)(u) in your fieldplot call.

And you also don't need to define W the same way, a second time. Since the original W was defined before the parameters were assigned values then calls to W or D(W) will pick up the latest parameter values assigned.

So your same processes can become more tidy and legible, eg, proteins_ac.mw

If you really want to visualize the latest parameter values in the formula then you can just invoke it to print it, eg.
   D(W)(Un);

If you want you can assign D(W) to some name, to avoid having D repeat effort. Eg,
   dW := D(W);
and then use that. Eg. proteins_acc.mw

You could put all your procedures into a package module, and then LibraryTools:-Save that instead.

See also the relevant section of the Programming Guide.

Here are various ways in which one may leverage the collect command to attain the target form.

Download collect_I_ex.mw

You can simply pass the explicit option, to get the result in terms of explicit radicals directly from the solve command.

Download solve_sys_expl.mw

In other words, you don't need to use two different commands to get the roots in radical form.

For a polynomial system involving more than one variable then results involving roots/radicals are not expressed in explicit form by default (even for this low degree).

See also the subsection Controlling the Form of Solutions in the Help page for topic solve,details . (Unfortunately, the difference between a single polynomial in one variable and system is not made entirely clear there, in this respect.)